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Abstract 

We formulate and study computationally the fluctuating compressible Navier-Stokes equations 
for reactive multi-species fluid mixtures. We contrast two different expressions for the covariance 
of the stochastic chemical production rate in the Langevin formulation of stochastic chemistry, 
and compare both of them to predictions of the chemical Master Equation for homogeneous well- 
mixed systems close to and far from thermodynamic equilibrium. We develop a numerical scheme 
for inhomogeneous reactive flows, based on our previous methods for non-reactive mixtures [K. 
Balakrishnan, A. L. Garcia, A. Donev and J. B. Bell, Phys. Rev. E 89:013017, 2014]. We study the 
suppression of non-equilibrium long-ranged correlations of concentration fluctuations by chemical 
reactions, as well as the enhancement of pattern formation by spontaneous fluctuations. Good 
agreement with available theory demonstrates that the formulation is robust and a useful tool in 
the study of fluctuations in reactive multi-species fluids. At the same time, several problems with 
Langevin formulations of stochastic chemistry are identified, suggesting that future work should 
examine combining Langevin and Master Equation descriptions of hydrodynamic and chemical 
fluctuations. 

PACS numbers: 05.40.-a,47.11.-j,47.10.ad, 47.70.Fw 

Keywords: Fluctuating hydrodynamics, Fluctuating Navier-Stokes equations, Multi-species, Thermal fluc¬ 
tuations 
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I. INTRODUCTION 


Chemical reactions are of central importance in both natural and industrial processes 
spanning the range of length scales from the microscopic, through the mesoscopic, and up 
to macroscopic scales. It is the rule, rather than the exception, that chemical reactions 
are strongly coupled to hydrodynamic transport processes, such as advection, diffusion, 
and thermal conduction. Prominent examples include diffusion-limited aggregation, pattern 
and chemical wave formation in reactive solutions, reaction-driven convective instabilities, 
heterogeneous catalysis, combustion, complex biological processes, and others. Even in a 
homogeneous system with only slightly exothermic reactions, the chemistry is coupled to 
the hydrodynamics, leading to non trivial effects such as fluctuation induced transitions [1]. 

Fluctuations affect reactive systems in multiple ways. In stochastic biochemical sys¬ 
tems, such as reactions inside the cytoplasm, or in catalytic processes, some of the reacting 
molecules are present in very small numbers and therefore discrete stochastic models are 
necessary to describe the system. In diffusion-limited reactive systems, such as simple co¬ 
agulation 2 A —> A 2 or annihilation A + B —)■ C, spatial fluctuations in the concentration of 
the reactants grow as the reaction progresses and must be accounted for to accurately model 
the correct macroscopic behavior. [2, 3] In unstable systems, such as diffusion-driven Turing 
instabilities [4-8], detonation [1], or buoyancy-driven convective instabilities [9], fluctuations 
are responsible for initiating the instability and may profoundly affect its subsequent tem¬ 
poral and spatial development. In systems with a marginally-stable manifold, fluctuations 
lead to a drift along this manifold that cannot be described by the traditional law of mass 
action, and has been suggested as being an important mechanism in the emergence of life 
[ 10 - 12 ], 

Much of the work on modeling stochastic chemistry has been for homogeneous, “well- 
mixed” systems, such as continuously stirred tank reactors (CSTRs), but there is increasing 
interest in spatial models [13]. When hydrodynamic transport is included, the focus has 
almost exclusively been on species diffusion, and there is a large body of literature on 
stochastic reaction-diffusion models. A Master Equation approach, notably, the Chemi¬ 
cal Master Equation (CME), is widely accepted for modelling well-mixed systems. The 
Reaction-Diffusion Master Equation (RDME) extends this type of approach to spatially- 
varying systems [14-16]. In the RDME, the system is subdivided into reactive subvolumes 
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(cells) and diffusion is modeled as a discrete random walk by particles moving between cells, 
while reactions are modeled using local CMEs [17, 18]. A large number of efficient and elab¬ 
orate event-driven kinetic Monte Carlo algorithms for solving the CME and RDME, exactly 
or approximately, have been developed with many tracing their origins to the Stochastic 
Simulation Algorithm (SSA) of Gillespie [19, 20]. The issue of convergence as the RDME 
grid is refined is delicate. [21-23] Variants of the RDME have been proposed that improve or 
eliminate the sensitivity of the results to the grid resolution, such as the convergent RDME 
(CRDME) of Isaacson [23] in which reactions can happen between molecules in neighboring 
cells as well. Particle-based spatial methods for stochastic chemistry include reactive Brow¬ 
nian dynamics [24, 25], Green’s Function Reaction Dynamics [26, 27], first-passage kinetic 
Monte Carlo [28-31], the small-voxel tracking algorithm [32], and others. 

In large part, the development of stochastic reaction-diffusion models has been divorced 
from work in the fluid dynamics community. Full hydrodynamic transport including ad- 
vection, sound waves, viscous stress, thermal conduction, etc., as well as nonequilibrium 
thermodynamics and chemistry, are fairly common in the reacting flow community. See, 
for example, textbooks by Kuo [33] and Law [34]. However work in this area focuses on 
macroscopic modeling; spontaneous thermal fluctuations, either chemical or hydrodynamic, 
are typically not considered. 

Within the field of nonequilibrium thermodynamics, the fluctuation-dissipation theorem 
provides the connection between hydrodynamic transport and spontaneous fluctuations. In 
particular, as an extension of conventional hydrodynamic theory, fluctuating hydrodynam¬ 
ics incorporates mesoscopic fluctuations in a fluid by adding stochastic flux terms to the 
deterministic fluid equations [35]. These noise terms are white in space and time and are 
formulated using fluctuation-dissipation relations to yield equilibrium covariances of the 
fluctuations in agreement with equilibrium statistical mechanics. Linearized fluctuating hy¬ 
drodynamics was first introduced by Landau and Lifshitz [36] and has since been used to 
study simple and binary fluid systems in and out of equilibrium [35]. 

A number of numerical algorithms for solving the equations of fluctuating hydrodynam¬ 
ics have also been developed [37-44]. These algorithms draw from a wealth of deterministic 
computational fluid dynamics (CFD) techniques and handle transport such as diffusion in 
a much more sophisticated fashion than random hopping between cells. For example, they 
include effects such as cross-diffusion, barodiffusion, thermodiffusion (i.e., Soret effect) as 
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well as advection by fluid motion. Furthermore, semi-implicit temporal discretizations and 
higher-order spatial discretizations can be used, even as fluctuations due to the discrete 
nature of the fluid are accounted for. In this spirit, Koh and Blackwell [45] propose using 
a more traditional gradient-driven diffusive flux formulation consistent with CFD practice 
within an CME-based description; however, their treatment of fluctuations is rather ad hoc 
and not consistent with the formulation of stochastic mass fluxes in fluctuating hydrody¬ 
namics. Very recently a spatial chemical Langevin formulation (SOLE) was proposed [46] in 
which the chemical Langevin approximation [47] is applied to the RDME treating diffusive 
hops as another reaction in a very large reaction network. While this leads to a formulation 
similar to fluctuating hydrodynamics it has several shortcomings, notably, it does not allow 
one to treat diffusion using advanced CFD algorithms. 

Investigations utilizing fluctuating hydrodynamics have revealed the crucial importance 
of hydrodynamic fluctuations in transport mechanisms, especially mass and heat diffusion. 
Notably, it is now well-known that all nonequilibrium diffusive mixing processes are accom¬ 
panied by long-range correlations of fluctuations. In certain scenarios these nonequilibrium 
fluctuations grow in physical extent well beyond molecular scales with magnitudes far greater 
than those of equilibrium fluctuations. These so-called “giant fluctuations” are observed in 
laboratory experiments [48-50], and arise because of the coupling between thermal veloc¬ 
ity fluctuations and concentration or temperature fluctuations. In fact, it has recently been 
shown using nonlinear fluctuating hydrodynamics that mass diffusion in liquids is dominated 
by advection by thermal velocity fluctuations [51]. Therefore, modeling diffusion using collec¬ 
tions of independent random walkers, as done in the RDME, is fundamentally inappropriate 
for describing the nature of hydrodynamic fluctuations at microscopic and mesoscopic scales; 
instead, hydrodynamic coupling (correlations) between the diffusing particles must be taken 
into account [52], Including fluctuations within the continuum description has also been 
shown to be important in particle-continuum hybrids [53], and should also benefit hybrid 
models for reaction-diffusion systems [54], 

In the hydrodynamic equations chemical reactions may be treated as a white noise source 
term [16, 55], in a fashion analogous to the stochastic transport fluxes. The study of fluc¬ 
tuating hydrodynamic models that include chemical reactions is relatively recent and there 
are few computational studies in the literature. Stochastic reaction-diffusion equations are 
considered by Atzberger in [56], but only within the reaction-diffusion framework and not 
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accounting for fluctuations in the chemical production rates. A thorough discussion of 
stochastic formulations of chemical reactions within the framework of statistical mechanics 
can be found in the monograph by Keizer [57]; Keizer does not, however, consider hy¬ 
drodynamic transport in spatially-extended systems in depth. In a sequence of important 
papers [58-61], chemical reactions have been incorporated in a nonlinear nonequilibrium 
thermodynamic formalism, making it possible to combine realistic nonlinear deterministic 
models based on the traditional law of mass action (LMA) kinetics with fluctuating hydro¬ 
dynamics. When considering fluctuations, however, a linearized approximation was used 
by the authors, limiting the range of applicability to modeling small Gaussian fluctuations 
around a macroscopic state that evolves in a manner unaffected by the fluctuations. A 
more phenomenological approach was followed to fit the LMA into the nonequilibrium ther¬ 
modynamics GENERIC formalism by Grmcla and Ottinger [62], but fluctuations were not 
considered. Here, we formulate a complete set of fluctuating hydrodynamic equations for a 
reactive multispecies mixture of ideal gases. We account for mass, momentum and energy 
transport, and chemical reactions, and consider a nonlinear formalism for describing the 
thermal fluctuations. 

Hydrodynamics is a macroscopic coarse-grained description, and fluctuating hydrody¬ 
namics is a mesoscopic coarse-grained description. As such, both descriptions rely on the 
approximation that the length and time scales under consideration are much larger than 
molecular, i.e, that each coarse-grained degree of freedom involves an average over many 
molecules. In fact, although formally written as a continuum model, fluctuating hydrody¬ 
namics is, in truth, a discrete model that only makes sense when seen as a coarse-grained 
description for the evolution of a collection of spatially-discrete hydrodynamic variables in¬ 
volving averages over many nearby molecules [63, 64], The fact that many molecules are 
involved in the reactions allows for a Langevin-like continuum description (i.e., diffusion 
processes) of the fluctuations instead of discrete models such as master equations (i.e., jump 
processes). The accuracy of Langevin formulations for chemically reacting systems has long 
been a topic of debate [65, 66]. In this work, we take the first step in combining realistic fluid 
dynamics with a stochastic chemical description and adopt a Langevin approach to describ¬ 
ing fluctuations. In future work, we will explore combining Langevin and ME approaches 
together, thus further bridging the apparent gap between the two. 

Here, we first formulate the fluctuating reactive Navier-Stokes-Fourier equations, discuss 
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their physical validity, and develop numerical methods for solving the stochastic partial 
differential equations. The methodology is a direct extension of our previous work on fluc¬ 
tuating hydrodynamics for non-reactive multispecies gas mixtures [42] to include a Langevin 
model of chemical reactions. We consider two distinct Langevin models, which are identical 
when very close to chemical equilibrium but differ far from thermodynamic equilibrium. 
The first model, which we term the Log-Mean equation (LME), is based on the GENERIC 
formulation of Grmela and Ottinger [62], but can be traced to older work on the subject 
as well [57, 67, 68]. The second Langevin model is the more familiar Chemical Langevin 
Equation [47, 57, 69]. 

The resulting algorithms are used to assess the importance of thermal fluctuations in 
several simple but relevant examples. The first example is a simple dimerization reaction, 
which has been studied theoretically in prior work by others [59-61]. Our second example is 
the Baras-Pearson-Mansour (BPM) reaction network [70, 71], which exhibits a rich behavior 
ranging from bistability to limit cycles. We study these examples in both well-mixed small- 
scale systems, comparing with the Chemical Master Equation, and in spatially-extended 
systems, comparing with fluctuating hydrodynamic theory and previous numerical work. 
For the latter, rather than imposing the non-equilibrium constraint by fixing concentrations 
in the bulk, the constraints are applied as boundary conditions, thus maintaining strict 
consistency with equilibrium thermodynamics, including microscopic reversibility (detailed 
balance), in all of the models we study. These examples illustrate how thermal fluctuations 
drive giant concentration fluctuations and how they affect the rate of pattern formation in 
an inhomogeneous system. 


II. THEORY 

In this section, we summarize the mathematical formulation of the complete fluctuating 
Navier-Stokes (FNS) equations for compressible reactive multispecies fluid mixtures. The 
details for non-reactive fluid mixtures are presented in [42]; here we focus on the chemistry. 
The formulation is first presented in its general form; the specific case of reactions in ideal 
gas mixtures is treated in Section IIC. 

The species density, momentum and energy equations of hydrodynamics for a mixture of 
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iVs species are given by 


d 

— ( p s ) + V • (p s v) + V ■ T s = m s n s , (s = 1,... Ns) (1) 

^ (pv) + V • [pvv r + pi] + V • n = pg, (2) 

^ (pE) + V ■ [(pE + p)v] + V ■ [Q + n • v] = pv ■ g, (3) 

where p s , m s , and are the mass density, molecular mass, and number density production 
rate for species s. The variables v, p, and E denote, respectively, fluid velocity, pressure, 
and specific total energy for the mixture. The total density is p = Ps-, g is gravitational 
acceleration, and superscript T denotes transpose. 

We consider a system with 7V R elementary reactions with reaction r written in the form, 

N S N S 

£H,: Y1 "J®** 

s= 1 

Here are the chemical symbols and are the molecule numbers for the forward and 
reverse reaction r. The stoichiometric coefficients are u sr — v~ r — uf r and mass conservation 
requires that 'Yh s u sr m r = 0. [57] For simplicity of notation, when there is no ambiguity we 
omit the range of the sums and write Y2 S ^ or sums over species, and write f° r sums 
over all reactions. Note that chemistry does not appear explicitly in the energy equation (3) 
since the species heat of formation is included in the specific total energy. 

Transport properties are given in terms of the species diffusion flux, viscous tensor, 
II. and heat flux, Q. Mass conservation requires that the species diffusion flux and the 
production rate due to chemical reactions satisfy the constraints, 

^ T s = 0 and m s Q s = 0 (4) 

S S 

so that summing the species equations gives the continuity equation, 

+ V ■ (pv) = 0. (5) 

The detailed form of the transport terms is summarized in Appendix A, see [42] for details. 
It is important to note that we neglect any possible effect of the chemical reactions on the 
transport coefficients of the mixture. 

We write the chemical production rate as the sum of a deterministic and a stochastic 
contribution, + h2 s , with the stochastic rate going to zero in the deterministic 



limit. [72, 73] To formulate these production rates we define the dimensionless chemical 
affinity as, 


A 


r 


-—— v sr m s u s 
k B T ^ r 


"srf 1 * ~ "srfrsr 


( 6 ) 


where fi s is the specific chemical potential (i.e., per unit mass) and ji s = m s ii s /k B T is the 
dimensionless chemical potential per particle; T and k B are temperature and Boltzmann’s 
constant, respectively. Summing over reactions gives the deterministic production rate for 
species s as [62] 

n» = Y' v sr f A r (7) 

■' r r k B l 

r 


where 



( 8 ) 


and T r is a time scale characterizing the reaction rate. This form of the deterministic equa¬ 
tions, while at first sight appearing different from the more familiar law of mass action 
(LMA), is fully consistent with it. The production rate, as given by (7) and (8), is also 
consistent with nonequilibrium thermodynamics [59]; this way of expressing the production 
rate in terms of a thermodynamic driving force (difference of exponentials of chemical po¬ 
tentials) can be seen as a generalization of the LMA to non-ideal systems, as elaborated in 
Section IIC. 

For a binary mixture undergoing a dimerization reaction, the deterministic part of the 
complete set of hydrodynamic equations including chemical reactions has been fit into a non¬ 
linear nonequilibrium thermodynamics formalism in Ref. [59] by introducing an additional 
reaction coordinate, as inspired by earlier work of Pagonabarraga et al. [58]. This extends 
earlier considerations of dimerization reactions in a strictly linear fluctuating chemistry for¬ 
malism [74], In the limit of high reaction barrier the equations written in [59] are equivalent 
to the ones we employ here even though our notation is different. However, fluctuating 
contributions in [59] are only considered in a linearized approximation, severely limiting 
the range of applicability to describing small Gaussian fluctuations around a deterministic 
average flow. 

In next two sub-sections we develop two nonlinear forms for the stochastic contribution 
to the reactive production rates, one coming from irreversible thermodynamics cast in the 
GENERIC formalism [62], and the other being a generalization of the more familiar form 
associated with the chemical Langevin equation (CLE) [47, 57, 69]. 
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A. The Log-Mean Equation 


Grmela and Ottinger [62] cast the phenomenological LMA (7) in the GENERIC formalism 
and obtain a nonlinear form for the dissipative matrix, under the assumption of a quadratic 
dissipative potential. Note that the entropy production rate can uniquely be written as 
a quadratic function of the thermodynamic driving force only for a single reaction; the 
resulting peculiar form of the mobility (dissipative) matrix (see Eq. (113) in [62]) involving 
a logarithmic mean has recently been justified from a model reduction perspective [75]. Here 
we assume that there is no cross-coupling between different reactions and thus associate 
an independent stochastic production rate with each reaction. Coupling between distinct 
reactions has been considered within a nonequilibrium thermodynamic framework only in 
some very specific cases [76, 77] and a general formulation requires more detailed knowledge 
about the coupling mechanism than is available in practice. 

Following the general principles for including fluctuations 1 in the GENERIC formalism 
[78], it is straightforward to write a Gaussian stochastic production rate assuming indepen¬ 
dence among the different reaction channels, 

^ Gr V^r M O Z? (9) 

r 

where 

pLM _ P -Ar 

TrksT Ar 

and where Zp 1 are independent white-noise random scalar fields with covariance 

(Zr( r, t)Z$( r', t')) = 8 r y S( r - r') S(t - t'), 

with each Z^ driving the stochastic production rate of a single chemical reaction r. We 
refer to this formulation for the stochastic chemistry as the “log-mean” form; the reasoning 
behind this name will become evident when presented in Section IIC for ideal mixtures. Note 
that (9) uses the kinetic or Klimontovich interpretation [79, 80] of the stochastic integral, 
formally denoted as a kinetic stochastic product with a o symbol in (9). The variance of the 
stochastic forcing XA M can be seen to be positive because A and A always have the same 
sign [62], Note that ^2 s m s tt s = = 0, as required by mass conservation. 

1 Fluctuations are not considered by Grmela and Ottinger [62], however, the “square root” of the mobility 
matrix written in Eq. (113) in [62] is straightforward to write, and this leads to the log-mean equation 
considered here. 
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For the purposes of exposition it is useful to consider a homogeneous “well-mixed” sys¬ 
tem of volume AV, which will correspond to a single hydrodynamic cell after spatial dis¬ 
cretization of (1). The dynamics of the (intensive) number density n s = N S /AV, where 
N s = AVp s /rri s is the number of molecules of species s in the cell, is given by, 



which can also be written in Ito form as, 



where W,^(t) are independent scalar white noise processes with covariance (W, sl (t)W^(f / )) = 
S r yS(t — t'). We call this system of stochastic ordinary differential equations (SODEs) the 
“log-mean” equation (LME). 

The derivative in the last stochastic drift term in (12) is the directional derivative of 
along the reaction coordinate. Unlike the more familiar Ito or Stratonovich interpretations 


of the noise, the kinetic form of the noise ensures that the corresponding Fokker-Planck 
equation has the traditional form [81], 


dP(n, t ) 
dt 


^ ^ dn s { UsrV ^ ArP + AV^ ‘ 

r s ! L s' 


This ensures that the LME is in detailed balance with respect to the Einstein distribution 
~ S/(ksT) for a closed system at thermodynamic equilibrium, where S is the total entropy of 
the system [78]. We note that it is not possible to obtain the LME from the chemical master 
equation (CME) with a systematic procedure; one must invoke some guiding principles about 
the structure of coarse-grained Fokker-Planck equations to “derive” this form of the noise 
[67, 68, 81]. 


B. The Chemical Langevin Equation 

Since both A and A are equal to zero at chemical equilibrium, near chemical equilibrium 
we can linearize (10) to first order in the affinity A, and approximate the amplitude of the 
stochastic production rate in terms of a sum over each forward and reverse reaction, that is, 




11 



Since this sum of products of exponentials is evidently positive, we can potentially use it 
even far from chemical equilibrium, and write the stochastic production rate as, 


= E ?y - v 2 ^ c + 


where 2 


v cl _ y 
Ur ± ~ 2 


P 


+ 


n 


z n_ 


^I'srft'S 


(14) 


(15) 


jj-ksT y 

Here Z^ + are independent white-noise scalar random fields that give the stochastic contribu¬ 
tion from the forward reaction, while Z\}_ correspond to the reverse reactions; the forward 
and reverse reactions are taken to be independent. In the next section the production rate 
factors, XA M and V^ L , are further simplified for the case of ideal gas mixtures. 

The form (14) for the amplitude of the stochastic production rate is found in most work 
on the subject [47, 59-61, 82], For example, though not written in this form, eqn. (8f) in 
Ref. [61] contains a sum of two exponential terms and is equivalent to (13) for the specific 
reaction considered there 3 . For a well-mixed homogeneous system of volume AH, the 
number densities of molecules of the different species follow the system of SODEs, 

dn s (t) 


dt 


£■ 


= > v, 


P 

TrksT 


A r + 


27b CL 
AH r 


(16) 


The stochastic equation (16) is commonly referred to as the chemical Lange vin equation 
(CLE) following Gillespie [47], and can be obtained from the CME by an uncontrolled 
truncation of the Kramers-Moyal expansion to second order. It is traditional to assume an 
Ito interpretation of the noise in the CLE, even though no precise justification for this can be 
made within the accuracy to which the CLE approximates the CME [69]. Mathematically, 
the nonlinear CLE contains similar information to the central limit theorem (i.e., linearized 
fluctuating hydrodynamics) corresponding to the CME in the limit of weak noise (large 
number of reactant molecules). 

As seen from (13), the two stochastic differential equations for the number densities, the 
LME using the kinetic noise (12) and the CLE using the Ito noise (16), are equivalent near 
chemical equilibrium. They are, however, different far from chemical equilibrium, as we 
2 Note that the two terms in (14) can be combined together to give a single white-noise process with 


amplitude I?, CL = this leaves the covariance of the stochastic forcing unchanged. 

3 As the authors note in a later publication [61], the sign in Eq. (88) in Ref. [59] is wrong and should 


be a plus rather than a minus. 


With this change that equation is equivalent to (14) evaluated at the 


deterministic solution, in the spirit of linearized fluctuating hydrodynamics. 
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illustrate in more detail in Section IV A. Notably, the forward and reverse reactions are 
treated together in the LME, consistent with the fact that, due to microscopic reversibility, 
there is only one independent rate coefficient for each reaction. [57] The ratio of the forward 
and reverse reaction rates is related to the equilibrium reaction constant, which is a thermo¬ 
dynamic and not a kinetic quantity. In fact, the LME is closely-related to the notion of the 
existence of a state of thermodynamic equilibrium in which each pair of forward and reverse 
reactions are in detailed balance with each other; one cannot write an LME for a system 
with irreversible reactions, which fundamentally violate detailed balance. 

By contrast, the forward and reverse reactions are treated completely independently in the 
CLE and there is no difficulty in writing a CLE for a system with irreversible reactions. The 
CLE is evidently inconsistent with the notion of detailed balance and is, in fact, inconsistent 
with equilibrium thermodynamics. Although written in a different form, Keizer’s (4.8.37) is 
the CLE, and Keizer’s (4.8.36) is the LME [57]; Keizer argues that the CLE is the correct 
equation and concludes: “Although the theoretical description of nonequilibrium ensembles 
would be greatly simplified if the phenomenological choice [LME] were correct, this appears 
not to be the case.” We will compare and contrast these two equations on some specific 
examples in Section IV A. 

C. The Law of Mass Action and Ideal Gas Mixtures 

In the formulation of hydrodynamic transport one normally works with the specific chem¬ 
ical potential, which has the general form, [83] 

ta s (p,T,X) = (\nX s + lnj s ) +n° s (p,T), 
m s 

where p° is the chemical potential at a reference state, X s — N s / N s > is the mole fraction, 
and 7 s (p, T, W) is the activity coefficient of species s. For chemistry it is more convenient 
to work with a dimensionless chemical potential per particle, 

hs ^7 ffi(Ais7s) + 

where jd° s = ( m s p,° s )/(ksT ). Note that X^g is the activity (i.e., effective concentration) and 
for an ideal mixture, r j s = 1. [84] This gives 

exp (vf r iis) = exp (vf r p,°) (X sls ) v ± r , 
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which leads to a generalized law-of-mass action (LMA) of the form 


ft, = 


r r k B T 


r \ s' 


{Xs'1 s 


fl A 


where nf(T,p,X) are the more familiar forward/reverse reaction rates (per unit time and 
per unit volume). Since there is only one independent timescale parameter, r r , the forward 
and reverse rates are not independent and the LMA gives the ratio to be the equilibrium 


= exp 


-£ 


constant, 

rv vn _ K r _ EL {X 8 lsY kr _ f y- 

Kr Ln.fe) -J eq V 3 J 

which is a purely thermodynamic quantity (closely related to the dimensionless reference 
Gibbs energy for the reaction at a unit reference pressure) that can be calculated from pure 
component data [82, 85]. Note that in chemistry texts the equilibrium constant is typically 
defined in terms of concentrations rather than activities as we have done here. 

For ideal gas mixtures we can further simplify the generalized LMA (17) to the more 
familiar form using number densities instead of mole fractions. From classical statistical 
mechanics, for an ideal gas mixture we can write 4 

p s = In n s + In > ( 19 ) 


fi s = In n s + In 


where A s = h/y/2jrm s kBT is the thermal wavelength of a structure-less particle and j(T) is 
the partition function for the internal degrees of freedom. [86] In general j(T) is a compli¬ 
cated function depending on the quantized energy levels of a molecule but in the classical 
approximation j(T) = (T/T 0 )^ z where z is the number of classical internal degrees of free¬ 
dom and T a is a reference temperature. 

For ideal gas mixtures the chemical production rate (17) can be written in the familiar 
power-law form, 


n, = 


V sr ( K IK- ~ k r I[n/ r ) ’ 


where kf are the forward/reverse reaction rates for the LMA formulated in terms of number 
density instead of activity. For uni-molecular reactions (e.g., QTTi —>■...) the “decay time” for 
a particle is usually assumed to be constant, in which case the corresponding reaction rate 

4 This is in the classical regime, that is, when the molecules’ mean spacing is much larger than their de 


Broglie wavelength. 
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(e.g., /c+) is a constant. For bi-molecular reactions (e.g., + WI 2 the production 

rate is usually assumed to be proportional to the collision frequency times an Arrhenius 
factor, in which case the corresponding reaction rate is only a function of temperature. [87] 
For the stochastic production rate in an ideal gas mixture, using, 


A r — y ] v sr ji s y ^ 


V sr Hs 


ln eX P (Eg V trfrs) 

eX P (Eg ’AvIA) 


In k nx- 

KU,nf r 


gives 


ls r 



In 

1 

2 


k n,< a+r -k n, n f r 
(k+ el n s tr ) - in (k el < 

KYl n f r +KYl n f r = 

s s 




( 21 ) 

( 22 ) 


where logmean and arthmean are the logarithmic and arithmetic mean, respectively. Note 
that T>\ M is zero if either reaction rate equals zero while "D,[? L is non-zero if cither reaction rate 
is non-zero. The logarithmic mean form of the noise in the ideal mixture case has appeared, 
phenomenologically, in several early papers [67, 68] and in more recent work [75, 88, 89] by 
other authors. 


III. NUMERICAL SCHEME 

The numerical integration of (l)-(3) is based on a method of lines approach in which we 
discretize the equations in space and then use an SODE integration algorithm to advance the 
solution using the basic overall approach described in [42], The spatial discretization uses a 
finite volume representation with cell volume AD, where U” ■ k denotes the average value of 
U = ( p s , pv, pE) in cell-(i, j, k ) at time step n. To ensure that the algorithm satisfies discrete 
fluctuation-dissipation balance, the spatial discretizations for the hydrodynamic fluxes are 
done using centered discretizations; see [40] and [42] for details. 

Discretization of the system in space results in a system of SODEs driven by a collection 
of independent white-noise processes W (t) that represent a spatial discretization of the 
random Gaussian fields 2(r,t) used to construct the noise. After temporal discretization 
these white noise processes are represented by a collection of i.i.d. standard normal variates 
Z, which can be thought of as a spatio-temporal discretization of Z\ the discretization is 
reflected in the presence of a prefactor 1 / \/ AV At in the expressions for given below [90]. 
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For temporal integration we use the low-storage third-order Runge-Kutta (RK3) scheme 
previously used to solve the single and two-component FNS equations [40], using the weight¬ 
ing of the stochastic forcing proposed by Delong et al. [90]. With this choice of weights, 
the temporal integration is weakly second-order accurate for additive noise (e.g., the lin¬ 
earized equations of fluctuating hydrodynamics [91]). As discussed at length in Ref. [91] 
the hydrodynamic stochastic fluxes should be considered as additive noise in a linearized 
approximation. 

The implementation of the methodology supports three boundary conditions in addition 
to periodicity. The first is a specular, adiabatic wall which is impermeable to momentum, 
mass or heat transport (i.e., all fluxes are zero at the wall). A second type of boundary 
condition is a no slip, reservoir wall at which the normal velocity vanishes (i.e., the total 
mass flux at the boundary vanishes) and the other velocity components, mole fractions and 
temperature satisfy inhomogeneous Dirichlet boundary conditions; this mimics a permeable 
membrane connected to a reservoir on the other side of the boundary. The third boundary 
condition is a variant of the no slip condition for which the wall is impermeable to mass 
but conducts heat. When a Dirichlet condition is specified for a given quantity, the corre¬ 
sponding diffusive flux is computed as a difference of the cell-center value and the value on 
the boundary. In such cases the corresponding stochastic flux is multiplied by y/2 to ensure 
discrete fluctuation-dissipation balance, as explained in detail [43, 92], 

Because the noise arising from the chemical reactions is multiplicative special care must be 
taken to capture the correct stochastic drift terms arising from the kinetic interpretation of 
the noise in the LME. [93, 94] We have chosen to write the equations in Ito form, which leads 
to an additional stochastic drift term in the LME (12). To integrate the Ito form in time, 
we evaluate the amplitude of the noise at the beginning of the time step and reuse the same 
random increments in all three stages of the RK3 scheme. The stochastic drift term arising 
in the LME is treated as a deterministic term but is also only evaluated at the beginning 
of the time step. The resulting scheme is only first-order weakly accurate. It is possible to 
construct second-order weak integrators by using the special one-dimensional nature (i.e., 
there is only a single reaction coordinate for each reaction even if there are many species 
involved) of the chemical noise [95]. However, in our simulations the time step is typically 
limited by stability considerations for advective and diffusive hydrodynamic processes, and 
therefore chemistry is accurately resolved even by a first-order scheme. Alternative temporal 
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integration strategies will be discussed in the Conclusions. 


The chemical Langevin form of the noise in (16) is discretized as, 


(ft 


CL\n 


H,j,k 




' 2 <»?%,* 


AVAt 


(K 




s = l,...,N s , 


(23) 


where , fc are zero-mean normal Gaussian variates generated independently in each 

cell at the beginning of each time step, ((Z^)” ■ lfc (Z^)” , ■, fc/ ) = 5i,ii8jj'5k,k'5n,n'5 r y ■ Note that 
the two terms in (14) have been combined into a single white-noise process with amplitude 
V ( r : Jj = V ^ + T>^_. As discussed above, the LME noise (9) leads to an Ito correction in 
( 12 ), 



The directional derivative of T)] M in the last term can be evaluated analytically, or, for 
simplicity of implementation, it can be efficiently approximated numerically using a finite 
difference along the reaction coordinate. 


IV. NUMERICAL RESULTS 

In this section we describe several test problems that demonstrate the capabilities of the 
numerical methodology. We consider two reaction systems, the first being simple dimeriza¬ 
tion, 

SHi : 2A <=> A 2 (24) 

where 5DT = (A,A 2 ), that is, species 1 is the monomer and species 2 is the dimer. In 
Section IV A 1 we investigate simple dimerization in a homogeneous system; in Section IV B 
we investigate the “giant fluctuation” phenomenon in the presence of dimerization for a 
system with an applied concentration gradient. 

The second model we consider is based on the Gray-Scott (GS) model [96, 97], which is 
known to exhibit a rich morphology of stationary and time-dependent patterns [98]. This 
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model, as formulated by Pearson [98], consists of the reactions, 


9ti 

: U + 2V -A 3V 


9*2 

: V ^ S 

(25) 

^3 

IT 


m 4 

: V^V f 


where the concentrations of the “feed species”, Uf and Vf, are held fixed and species S is 

inert. Since elementary reactions are 

rarely trimolecular in nature, 

we consider a variation 

of the GS model developed by Baras et al. [70, 71]. The Baras-Pearson-Mansour (BPM) 

model 5 is, 



Dir : 

U + W ^V + W 


9l 2 : 

2V^W + S 


91 3 : 

Zr> 

IT 

(26) 

91 4 : 

IT 


91. 5 : 

IT 



with 9Jt = (U,V,W, S,Uf,Vf). This model was developed as a more realistic variant of 
the Gray-Scott model suitable for particle simulations of dilute gases. 6 Note that the first 
and third reactions are irreversible in the original BPM model, which is not consistent with 
detailed balance. The BPM model has not been studied as extensively as the GS model but 
its dynamics are expected to be qualitatively similar. Note that the BPM model replaces the 
trimolecular reaction in the GS model with a pair of bimolecular reactions and introduces 
W as an intermediary species. In the standard BPM model [70] the number densities of 
S, Uf , and Vf are held fixed so, being an open system, total mass is not conserved and 
detailed balance is not satisfied. [57] In Section IV A 2 we investigate this standard BPM 
model in a homogeneous “well-mixed” system; in Section IV C we simulate a two dimensional 
domain with full hydrodynamic transport with species S, Uf, and Vf held fixed only at the 
boundaries. 

5 Also known as the OLR model. 

6 To implement the BPM model in a molecular simulation with binary collisions Baras et al. modify the 
model by making all the reactions bimolecular and by introducing an “auxiliary” species, A. 
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A. Homogeneous Systems 


We first consider homogeneous “well-mixed” 7 systems of volume AV with only chemistry 
(i.e., no hydrodynamics). In this section we compare the results obtained using the log-mean 
equation (LME) form, (21), and the chemical Langevin equation (CLE) form, (22), with 
results from CME simulations performed using the Stochastic Simulation Algorithm (SSA), 
also known as the Gillespie algorithm. [19] The chemical master equation (CME) is widely 
accepted as an accurate model for well-mixed chemical systems and SSA is a popular scheme 
for simulating the stochastic process described by the CME. [20] As we will see, the two 
forms for the Langevin noise have their advantages and disadvantages and both forms are 
only approximations of the CME with limited ranges of validity. 


1. Dimerization Reaction 


We start by considering the dimerization reaction (24) in a closed system for which the 
deterministic production rate for species 1 (monomers) is 

fli = — 2{k + n\ — /c _ n 2 ), 

and by mass conservation, for dimers fl 2 = — since 'm 2 = ‘2m 1 . The constraint of mass 
conservation can also be expressed as ri\ + 2n 2 = no, where no is the initial number density 
of A particles (in either monomer or dimer form). We may then write, 

fli = —2 k + n\ + k~(rio — n\) 


and limit our attention to the monomer species. For simplicity we take the ratio of the rate 
constants to be k + /k~ = 1/no so the equilibrium mass fraction (Yi) eq = ^ (i.e., Oi = 0 for 

Yi = !)■ 

The log-mean stochastic production rate for ri\ is 


O, = -2 


2V lm 

AV 




with 


'pLM _ 


k + n\ — k n 2 


In {k + n\) — In (/c _ n 2 ) 


7 The precise mathematical definition of a “well-mixed” system is delicate and will not be discussed here 
at length. Roughly speaking, it means that diffusion is sufficiently fast compared to reactions; see [73] for 
a theorem on the limit of infinite diffusion. 
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From the corresponding Fokker-Planck equation (FPE) (II A) one finds the LME is in de¬ 
tailed balance with respect to the equilibrium distribution 

fill = Z - 1 exp |»„AV ( lull - Vi) tv, - ll - r, In il')) - ( Y , (In (2) - 1) l, 

(27) 

where Z is a normalization constant. This Einstein distribution Pj q M (Y]) ~ exp(S (Y\) / ks) 
is in agreement with the the correct thermodynamic entropy S in the limit of Stirling’s 
approximation, as we demonstrate in Appendix B. 

For the chemical Langevin equation the stochastic production rate for ri\ is 

~ /9T)CL i 

fli = — 2 Y ^ yV n with V CL = - ( k + n\ + /c“n 2 ) . 

The equilibrium distribution can also be found from the stationary solution of the FPE 
corresponding to the CLE, which we do not write here for brevity. 8 We do note that, for 
this example, -P e q L (Pi) is quite close to a Gaussian. We further observe that, unlike the 
LME, no thermodynamic interpretation can be given to In P^. In fact, the tails of are 
quite different from those of P C q M and, being nearly Gaussian, the former includes unphysical 
values of the concentration (i.e., Y\ is not constrained between 0 and 1). 

Figure 1 shows numerical results for the equilibrium distribution of the number of 
monomers N { = riiAV = pY\AV/ni\. At thermal equilibrium the simulation results using 
the log-mean equation (LME) form for the noise are in excellent agreement with equilibrium 
statistical mechanics (see Appendix B) and with CME/SSA simulation results. Other work 
has also shown that, when detailed balance is obeyed, the LME correctly reproduces the 
equilibrium transition rates for rare jumps between stable minima in bistable systems [68]. 
On the other hand, the chemical Langevin equation (CLE) result has the noticeable flaw 
that, being a Gaussian, the distribution extends to unphysical negative values of N\. 

However, the LME does not compare favorably with the numerical solution of the CME 
for time-dependent situations, such as when a system is relaxing toward equilibrium. To 
illustrate this, we simulate an ensemble of systems prepared with an initial condition far from 
equilibrium, specifically with Y\ (t = 0) ~ 1, and measure the time-dependent probability 
distribution P(Tj,t) as the system relaxes toward chemical equilibrium ((Tj) eq = |). As 
expected, for the ensemble mean value of the number density hi(t) = (ni(t)}, we find close 

8 In the literature this equation is called the CFPE (chemical FPE), see (4) in [65]. 
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FIG. 1. Empirical histograms of the probability distribution P(N\) for the number of monomers Ah 
obtained from the log-mean equation (left) and the chemical Langevin equation (right), compared 
to results from the CME (obtained using SSA). The theoretical solution of the corresponding 
FPE is shown with a solid red line; the numerical results from SSA are the circles. The Einstein 
distributions employing the Gaussian approximation and the second-order Stirling approximation 
(B4) to the true entropy (Bl) are also shown. Note that (Ah) ~ 16, which is rather small and thus 
tests the limits of applicability of a Langevin (non-discrete or non-integer) description. 


agreement among LME, CLE, and CME results (not shown), even when fluctuations are 
quite large. However, if we consider the standard deviation of the number of monomers, 
the left panel in Fig. 2 clearly demonstrates that the CLE is in much better agreement with 
the CME (as shown by the SSA results) for describing relaxation toward equilibrium. Also 
shown on this graph is the theoretical solution for the standard deviation obtained by first 
linearizing the CLE 9 around the solution of the deterministic law of mass action (which 
is the law of large numbers corresponding to the CME [69]), and then writing a system of 
ODEs for the mean and variance of n\(t). Specifically, we have that dfi\(t)/dt = fA and, 
using Ito’s formula, we get the central limit theorem corresponding to the CME [69], 

: D cl (fii(t)), 


dt 

where C\(t) = ((ni(t) — h 1 (t)) 2 ). 


dCi (t) = ^ Mg) + 


dn i 


AV 


9 The linearized CLE is known as the second-order van Kampen expansion [99] in the physics literature 
and has the same formal order of accuracy as the CLE [69] but is much simpler to solve analytically due 
to its linearity. Analysis in Ref. [65] suggests that the nonlinear CLE may be more accurate than the 
linearized CLE for large noise but we are not aw^e of rigorous mathematical estimates. 










FIG. 2. Relaxation toward equilibrium in a closed cell initially containing 100 monomers and 4 
dimers and relaxing toward the equilibrium state of (Ni) ~ 54 monomers. Rates are k~ = 0.3 and 
k + = 2.78 • 10 , time step is At = 0.005. (Left) Time evolution of the standard deviation of the 

number of monomers; theory curves are obtained by solving (IVA1). (Right) Histogram of the 
probability density P(N\, t) for the number of monomers at an early time t = 0.05. 


The agreement between CLE and CME in the left panel of Fig. 2 is not surprising 
since it is well-known that the central limit theorem for the CME is a linear Langevin 
equation with the noise covariance given by the CLE form rather than the LME form. [69, 
100] The agreement between the CLE and SSA results is less impressive when we look 
more closely at the probability distribution during the relaxation toward equilibrium. The 
right panel of Fig. 2 shows histograms of the probability distribution for the monomers, 
P(Ni,t), at an early time in the relaxation. This distribution is not close to Gaussian 
for the CME/SSA solution and we see that, in this regard, the CLE result is no better 
than the LME result. In fact, for the probability distribution of the dimers (not shown) 
the CLE results have the unphysical feature of non-zero probability for negative values of 
dimer concentration. Of course, one can argue that the number of molecules in the system 
is too small for a Langevin approximation to apply. If the fluctuations are decreased, the 
probability distribution P(Ni,t) will become closer to Gaussian and then the CLE will 
provide a better description; note however that the tails of the distribution will always be 
incorrect for the CLE, even at thermodynamic equilibrium. 
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Species 

S 

Uf 

Vf 


n r (fixed) 

7.0 • 10 2 

5 • 10 11 

5.65 • 10 10 

- 

9*i 

^2 

9*3 

9*4 

5*5 

k + 

2•10” 3 

10 -3 

0.0200936 

0.28 

0.28 

k~ 

2•10~ 9 

2.198 • 10” 4 

2.009 • 10" 8 

2.8 • 10“ 10 

2.8 • 10“ 10 


TABLE I. Table of parameters for the BPM model in a well-mixed system (see Figure 3). The 
system volume AV = 4 and the time step At = 1.0 (reducing the time step size further did not 
change the results significantly). 


2. Bistable BPM Model 

As a less trivial homogeneous example, we consider the BPM model (26) for an open 
system held at a non-equilibrium steady state for which the probability distribution function 
is bimodal. [71] For the chemistry-only study in this section, the number of molecules of 
species 1, 2, and 3 (U, V, and W ) are allowed to vary while the number of molecules of all 
other species are fixed. The relevant parameters are given in Table I. Note that a similar 
system (with all reactions being bimolecular) was studied by Baras et ah, who found good 
agreement between SSA and molecular simulations using the Direct Simulation Monte Carlo 
(DSMC) algorithm [71]. Baras et al. also examined the accuracy of the CLE linearized 
around the solution of the deterministic equations, and, not surprisingly, found it to be a 
very poor approximation of the CME for the parameters they chose. Gillespie [47] suggests 
that “A repetition of the study of Baras and co-workers using the Langevin equation [CLE] 
instead of the [linearized CLE] should show the chemical Langevin equation in a fairer light.” 
This section presents such a study using both the CLE and the LME. The parameters were 
selected such that the number of particles is large (roughly O(10 2 — 10 3 ) for each species) 
but not so plentiful as to prevent the SSA simulation from accurately sampling the bimodal 
distribution in a reasonable amount of computation time. 

A phase-space picture of a typical trajectory is shown in the left panel of Fig. 3. The tra¬ 
jectory moves between two basins centered around the two stable deterministic steady states, 
which are labeled state A (corresponding to Ni ~ 1740, N 2 ~ 448, N 3 & 328 molecules) and 
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FIG. 3. Numerical results for the BPM modelin a well-mixed system (see Table I). (Left) A phase 
space picture of a typical trajectory from SSA as it visits the two basins around the deterministic 
steady states (labeled A and B in the figure). The insert shows the trajectory in the collective 
coordinate, x(t), based on assigning each point to either state A (red) or state B (green) based on 
the last state visited by the trajectory (see text). (Right) Histogram of the steady-state probability 
distribution P(x) comparing SSA, LME, and CLE simulation results. Measurements were skipped 
for the first 5 ■ 10 5 steps to relax the initial transient and then statistics were collected for 1 • 10 8 
steps. 


state B (corresponds to N± & 1224, N 2 « 936, N$ « 1424 molecules). Based on this picture, 
we chose to define a collective coordinate x(t) which is the projection of the state (ni, n 2 , n 3 ) 
onto the line connecting the two stable points (red line in the figure). This simple linear 
collective coordinate has the property that x = 0 at state A and x — 1 at state B ; note that 
x(t) is not bounded between zero and one. The insert in Fig. 3 (left panel) shows x(t ) for a 
typical trajectory as the system moves between the basins. 

The right panel in Fig. 3 shows the steady-state probability distribution P(x ) for the 
collective coordinate x, clearly illustrating its bimodal form; similar results are found for the 
probability distributions for ri \, n 2 , and 7i 3 . The results for the LME and CLE Langevin 
approximations are qualitatively similar to those from the CME/SSA but quantitatively 
different; the LME result is in better agreement with the CME for this specific example. 
To also examine the long-time dynamics of the well-mixed bistable BPM system, we assign 
each (discrete) point in time to either state A or state B (see insert in left panel of Fig. 3). 
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State A 

Mean 

Variance 

State B 

Mean 

Variance 

SSA 

1.13T0 5 

1.53-10 10 

SSA 

2 .86-10 5 

8.74-10 10 

LME 

3.23-10 5 

9.03-10 10 

LME 

8.45-10 5 

1.91-10 11 

CLE 

1.38-10 5 

1.70-10 10 

CLE 

5.01-10 5 

7.91-10 10 


TABLE II. Numerical results for the mean and variance of the time spent in state B before transiting 
to state A, and vice versa, for the three implementations of stochastic chemistry. 

The assignment is performed by defining two sets A = {x < 0.3} and B = {x > 0.6} and 
assigning each point in the trajectory to the last set that was visited. The distribution of 
waiting times spent in the two states before transitioning to another state is related to the 
transition rate, and the ratio of the average waiting times gives the ratio of the probabilities 
to be found in each of the two states. For large AV (weak noise), the transitions are rare 
events and the distribution of waiting times should be approximately exponential (recall 
that for an exponential distribution the variance is the square of the mean). Numerical 
results for the mean and variance of the time spent in state B before transiting to state A, 
and vice versa, are shown in Table II. 

Our results indicate that for the BPM model the CLE and LME provide a reasonably good 
qualitative description of the long-time dynamics and rare-event statistics for the parameters 
studied here. However, both approximations are in general uncontrolled and the quantitative 
match between the CME and either CLE or LME will not improve even if the cell volume 
increases and the fluctuations decrease in amplitude. In Ref. [68] it is observed that the 
LME correctly reproduces the very long-time dynamics (more precisely, the large deviation 
theory) of the CME for the bistable Schlogl model. This conclusion is, however, specific to 
this simple one-dimensional model because the system obeys detailed balance even though it 
is not in thermodynamic equilibrium; the BPM model studied here is not in detailed balance 
and there is no a priori reason to expect the LME to be more accurate than the CLE. The 
fact that the CLE is not able to describe rare events is well-known, see for example the 
discussion by Gillespie in [101] and after Eq. (9b) (which is the CLE) in [66], or recent 
numerical studies of noise-induced multistability [102], It can, in fact, easily be proven that 
this problem is shared by all diffusion process (SODE) approximations of the CME 10 , and 
10 While Langevin approximations cannot approximate atypical statistics of the CME, the CLE is likely 
appropriate for describing the typical behavior of the CME [65]. 
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fundamentally stems from the difference between the rare-event statistics of Gaussian and 
Poisson noise 11 . A promising alternative is to use tau-leaping to approximately integrate 
the CME in time [66] since it uses Poisson noise, and thus has the potential to correctly 
approximate the long-time behavior of the CME. This point is discussed further in the 
Conclusions. 

B. Giant Fluctuations 

We now consider a system for which concentration fluctuations are affected by both chem¬ 
istry and hydrodynamics in an interesting fashion. In the absence of chemistry a gradient of 
concentration induces a long-ranged correlation of concentration fluctuations [35, 104, 105]. 
These correlations are closely tied to the experimentally observed “giant fluctuation” phe¬ 
nomenon [48-50]. In an isothermal, nonreacting binary mixture the static structure factor 
for fluctuations in the mass fraction of the first species contains two contributions, 

S (k) = ((<W 1 )(5? 1 )*) = S eq + S neq , 

where “hat” denotes a Fourier component; the equilibrium part is 

= y(U)„ (1 - (U)eJ , (28) 

The non-equilibrium enhancement of the static structure factor due to a concentration gra¬ 
dient is S neq (k) ~ (Vli) 2 //c 4 , where the wavevector k is perpendicular to the imposed 
concentration gradient. 

The nature of these long-ranged correlations is modified in the presence of chemical reac¬ 
tions, as predicted by linearized fluctuating hydrodynamics [60, 61, 106]. Some preliminary 
numerical studies of fluctuations in the presence of chemistry have been done in Ref. [107] 
using an RDME-based description. However, these simulations are for a simpler isomeriza¬ 
tion A ^ B in one dimension and, furthermore, they are concerned with reaction-diffusion 
only and do not account for the hydrodynamic velocity fluctuations that are responsible for 
the giant concentration fluctuation phenomenon. 

We consider here the dimerization reaction (24) in a spatially inhomogeneous system. 
A rather detailed linearized fluctuating hydrodynamic theory for this example has been 
11 To see this contrast the quadratic Hamiltonian in (1.4) in [103], which applies to SODEs, and the expo¬ 
nential Hamiltonian in (1.6) in [103], which applies to master equations. 
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developed by Bedeaux et al. in [61], for a system in which a concentration gradient is 
imposed via a temperature gradient through the Soret effect. However, this analysis assumes 
a liquid mixture (large Schmidt and Lewis numbers) and thus does not apply to gas mixtures. 
Therefore, a simplified theoretical analysis of giant fluctuations in binary gas mixtures in 
the presence of an imposed constant concentration gradient and reactions is developed in 
Appendix C. 

Our simplified theory decouples the temperature equation and uses a concentration equa¬ 
tion (specifically the mass fraction of the first species) coupled to an incompressible fluctuat¬ 
ing velocity equation. For the case of a liquid mixture with very large Schmidt number, which 
is the case considered in [61], the calculation predicts that the nonequilibrium enhancement 
of the static structure factor of concentration fluctuations for the monomer species is (see 
eqn. (Cl)) 

(fc) = hTAW (1 + WT 1 • (29) 

where D is the diffusion coefficient and rj is the viscosity. The last term on the r.h.s. depends 
on the penetration depth d [108], 



We see that for large wavenumbers (k 3> 1 /d) the spectrum is ~ k ~ 4 , as in the absence of 
the chemical reaction. However for small wavenumbers (k <C 1/d) there is a transition to a 
k~ 2 spectrum. For gas mixtures, however, a more detailed model is required that takes into 
account the finite value of the Schmidt number. The result of this calculation is eqn. (C2), 
which predicts a further transition to a flat (constant) spectrum at small wavenumbers, with 
a finite S neq (k = 0) = (. ksT) (VYi) 2 /(9p(k~) 2 ). The calculation in Appendix C indicates 
that this effect is important even in liquids and the more refined theory ought to be used if 
quantitative agreement with experiments or simulations is sought. 

We performed a series of simulations to investigate these predictions using the full hy¬ 
drodynamic equations plus chemistry. It is important to note that even though we use the 
full nonlinear equations, nonlinearities in the fluctuations are negligible for the simulations 
reported here 12 . In fact, the noise is very weak (since the domain is quite thick in the z 
direction) and the numerical method is effectively performing a computational linearization 
of the fluctuating equations around the solution of the (nonlinear) deterministic equations 

12 The nonlinearity of the deterministic (macroscopic) equations is fully taken into account 
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[91]; this is simular to what is done analytically in Refs. [60, 61] but does not require any ap¬ 
proximations. In the small Gaussian noise regime the linearized CLE equation applies, and 
therefore in these simulations we use the CLE form for the stochastic chemical production 
rate in agreement with the theory in [60, 61]. Identical results (to within statistical errors) 
are obtained by using the LME form of the noise (not shown); this is not unexpected since 
the important noise here is the stochastic momentum tensor driving the velocity fluctua¬ 
tions; the stochastic mass flux and production rates only affect the reaction-diffusion part of 
the spectrum, which is much smaller than the nonequilibrium enhancement we study here. 

Here we assume that the traditional number-density based LMA (20) holds with constant 
rates k + and k~ . This requires that the time scale for the reaction is proportional to the 
number density, that is, r ~ p/ksT = n. From (18) and (19), for the dimerization of an 
ideal gas the ratio of these rates is, 

k+_ = (A 3 ) 2 ]2 = 2 3/2 A 3 32 

k- At tA) 2 ) 2 

which is a complicated function that depends on the form of the internal degrees of excitation. 
These details determine the number fraction (AR) eq or, equivalently, the mass fraction (Lj ) eq 
at chemical equilibrium; here we set the ratio of the forward and reverse rates to ensure 
(Yi)e q = 1/2, assuming that Ai, j\ and A are consistent with this choice. Since the reaction 
here changes the number density and thus the pressure, the reaction is strongly coupled 
to the momentum and energy transport equations. In order to minimize this coupling, we 
adjust the number of internal degrees of freedom of the dimer particles. Specifically, we set 
the heat capacities to c P: i = ^ks/mi (corresponding to three translational and zero internal 
degrees of freedom), c Pi2 = Sks/m-i (corresponding to three translational and five internal 
degrees of freedom). This choice ensures that c v of the mixture is independent of composition 
so that at constant pressure the reaction is isothermal. 

The fluid was taken to be a dilute binary mixture of hard-sphere gases, using kinetic 
theory formulae for the transport coefficients [39]. In CGS units, the species diameters are 
(j 1 = 2.58 • 10“ 8 and a 2 = 3.23 • 10 -8 , and m 1 = 6.64 • 1CT 23 . At equilibrium the density 
p eq = 1.78 • 1CT 3 , the temperature T eq = 300, and the concentration (yj) eq = 0.5. For these 
parameters, at the equilibrium conditions the mass diffusion coefficient is D = 0.2698 while 
the momentum diffusion coefficient (kinematic viscosity) is v = 0.2374. The ratio of the 
reverse and forward reaction rates is fixed at k~/k + = p/m 1 = 2.67985 • 10 19 . We vary the 



penetration depth d by changing the value of the reaction rates. 

The simulations used a 128 2 grid and time step size At = 2.5 • 10 -8 , grid spacing Ax = 
Ay = 10 -3 , and thickness in the z direction of Az = 10 -3 . The first 6 • 10 4 time steps 
were skipped and then statistics collected for 2 • 10 6 time steps. A steady concentration 
gradient was imposed by using Dirichlet boundary conditions at top and bottom boundaries 
(y — L and y = 0). Specifically, we take Y (y = 0 ,t) = 0.3 and Y\ (y = L y ,t) = 0.7, with 
temperature fixed at T = 300 and no-slip boundary conditions for the velocity. Periodic 
boundary conditions were used in the other direction. Concentration profiles for various 
values of penetration depth, d, are shown in Fig. 4. As expected, when the chemistry is 
slow (d: L) the concentration profile is nearly linear; when the chemistry is fast the 

concentration is nearly constant (at its chemical equilibrium value of (Yi) eq = 1/2) except 
near the boundaries. Note that we set the thermal diffusion ratio to zero (i.e., no Soret 
effect) so that the system is isothermal and the simple theory presented in Appendix C 
applies. 

For comparison, we also performed simulations in which we turn off all hydrodynamics 
except Fickian diffusion, giving us a reference reaction-diffusion structure factor S r d(k). For 
the case of a binary mixture [39] with k~/k + = p/m\ , the reaction-diffusion CLE reduces to 


£ (Y.) = V • (DVY0 + k~ (—2Y/ + (1 - Y,)) 


+ V • 



Y,(l - Y, 2 ) Z* 



(2Y] 2 + (1 - Y)) Z n , 


(30) 

(31) 


where n 0 = p\jrn\ + 2p 2 /'^2 is the total number density of A particles contained in both 
monomers and dimers, which is spatially constant in this reaction-diffusion approximation. It 
is important to note that the reaction-diffusion model (31) is thermodynamically inconsistent 
because it ignores the coupling of the chemistry to the energy and momentum transport. It 
is well-known that adding chemistry should not change the fluctuations at thermodynamic 
equilibrium [60], and this is indeed the case for the complete set of hydrodynamic equations 
that we study in this work. By contrast, for the reaction-diffusion (31) the equilibrium 
structure factor for (Yi) eq = 1/2 is given by 

( rd ) _ 3mi 8/9 + k 2 d 2 

eq 8p 1 + k 2 d 2 ’ 1 ’ 

which only approaches the thermodynamically correct answer (29) for kd 1. Note that 
the inconsistency between full hydrodynamics and reaction-diffusion is not evident in Eqs. 
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FIG. 4. Average concentration profile Y\ (y) at steady state for various values of penetration 
depth, d. Lines are from simulations using the full hydrodynamic equations; symbols are for 
reaction-diffusion only. 

(27a,28) in [60], because the authors of that work “neglect the dependence of the specific 
Gibbs energy difference on pressure.” This inconsistency is not of any importance in our 
study because we only use the reaction-diffusion simulations to obtain a baseline to subtract 
from the full hydrodynamic runs at large wavenumbers; at small wavenumbers the nonequi¬ 
librium enhancement is many orders of magnitude larger than the difference between (29) 
and (32). 

The reaction-diffusion runs are not limited by the Courant condition so we increased the 
time step to At = 2.5 • 1CT 7 ; a total of 6 • 10 4 steps were skipped initially and then statistics 
were collected for 3-10 6 steps. As seen in Fig. 4 the average concentration profiles are nearly 
the same whether the simulations used the full hydrodynamic equations or simply species 
diffusion. For the structure factor, however, we find that the reaction-diffusion simulations 
do not reproduce the giant fluctuation result (29), rather, they follow (32), which does not 
show a power-law enhancement, as seen in the top panel of Fig. 5. This is expected since 
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the giant fluctuation effect arises from the coupling of concentration fluctuations with the 
velocity fluctuations. 

Because the equilibrium structure factor (28) is derived assuming a uniform bulk state, 
which is not actually the case here, we define S neq (k ) = S(k) — S r d[k) as a measure of the 
“giant” nonequilibrium fluctuations coming from the coupling with the velocity equation. 
Results for S neq (k ) for several penetration depths are shown in the bottom panel of Fig. 5, 
comparing simulation results with the simple theory, eqn. (C2). Since chemistry should have 
minimal effects for large k according to the theory, we compute an effective concentration 
gradient by approximately matching the tail of the numerical result to the tail of the theory. 
We see that the theory correctly reproduces the qualitative trends, namely, that the giant 
fluctuations level off to a constant value at a wavenumber of order d~ 1 . However, except for 
the case of no reaction (d —>• oo), 13 the theory is not in quantitative agreement with the 
simulations. To confirm that the issue is not under-resolution of the penetration depth by the 
grid, we perform runs with a finer grid of 256 2 cells 14 , and we get the same result over the 
common range of wavenumbers, showing these runs are sufficiently resolved for the purpose of 
computing S[k). Note that in the plots the numerical wavenumber k x is corrected to account 
for discretization artifacts in the standard 5-point Laplacian, k 2 = sm 2 [k x Ax / 2) / {k x Ax / 2) 2 . 

The mismatch between theory and simulation is not so surprising since the theory is for a 
weak gradient applied to a system that is essentially near equilibrium; this is not true in this 
setup. The only way to get this isothermal system out of equilibrium is via the boundaries, 
so the system is actually far from chemical equilibrium near the boundaries and then goes to 
chemical equilibrium in the middle of the domain, but in the middle the gradient disappears. 
A new more sophisticated theory is required that linearizes not around a constant state but 
rather around a spatially-dependent state (this is automatically done in our code). Also, 
boundaries (i.e., confinement effects) may need to be included, especially for penetration 
depths comparable to the system size. 


13 Note that in this case the mismatch between the theory and simulations at very small k is due to the 
effect of boundaries [109]. 

14 For the more resolved runs Ax = Ay = 5 • 10~ 4 , At = 6.25 • 10~ 8 for the problem without hydrodynamics 
and At — 6.25 • 10~ 9 with hydrodynamics. 
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FIG. 5. (Top) The structure factor S r d(k) for the reaction diffusion model (31) in the presence of 
an applied gradient. For small d (fast reaction) the fact that the structure factor is not perfectly 
flat (constant) can be explained by (32) and comes from the thermodynamic inconsistency of the 
reaction-diffusion model. (Bottom) Non-equilibrium structure factors S neq (k) = S(k) — S r d(k ) for 
values of penetration depth d varying from oo (no reaction) to d = 10 3 , the width of one grid cell. 
Numerical results are shown with symbols, and the theoretical prediction (C2) is shown as a line 
with the same color. In the absence of a gradient is flat, S eq ~ 1.4 • 1CF 20 . 
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C. Pattern Formation 


Since the seminal work of Turing [110], pattern formation in deterministic reaction- 
diffusion systems has been investigated extensively, mostly in theoretical studies but also by 
laboratory experiments [111, 112]. The study of stochastic systems is more recent and, as 
described in the introduction, has primarily focused on models based on a reaction-diffusion 
master equation (RDME). Such a model was introduced by Mansour and Houard [113] as 
a practical numerical scheme for the study of correlations in spatially-distributed reactive 
systems. Subsequently RDME-based models have been used to study the influence of fluc¬ 
tuations on pattern formation for a variety of reaction-diffusion systems [4-7]. Recently, a 
spatial chemical Langevin formulation (SOLE) was proposed [46] and the RDME, SOLE and 
deterministic equations were compared for the development of patterns in the Gray-Scott 
model; it was observed that the SOLE is qualitatively similar to the RDME for the majority 
of examined sets of parameters, but not always. All of the RDME-based models usually use 
simplified descriptions of diffusion, but recently it has been observed that accounting for 
cross-diffusion effects (which are included in complete generality in our formulation) may 
lead to qualitatively-different behavior for Turing instabilities [114-117]. Particle simula¬ 
tions including full hydrodynamics have also been performed using the DSMC method [118] 
and molecular dynamics [8]; these are, however, limited to small systems in (quasi) one 
dimension because of the high computational cost of particle simulations. 

In our final example we consider pattern formation in the Baras-Pearson-Mansour (BPM) 
model (26) for a dilute gas mixture with full hydrodynamic transport. The system is initial¬ 
ized in a uniform constant “reference” state in which the number densities of the different 
species are as specified in Table III. These number densities and the reaction rates are set 
so that the reference state is similar to that investigated in [70]. Specifically, under the 
assumption that the number densities of the reservoir or “solvent” species S, Uf, and Vf are 
fixed, the deterministic dynamics of the three reactive species I/, V, and W starts close to 
the single unstable fixed point; the stable attractor of the dynamics is a limit cycle around 
this unstable point. Of course, when the number densities of the solvent species are not 
fixed the dynamics is six-dimensional and much more complex. Since we consider a time- 
dependent non-equilibrium scenario the chemical Langevin form of the noise (23) was used 
for the stochastic chemistry. 
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In the standard BPM model the solvent species (S, Uf, and Vf ) have hxed concentrations 
but in our hydrodynamic simulations they were hxed only at the boundaries. These three 
species are also made abundant to buffer them from having rapid variations in concentration 
(see Table III). In the limit of infinite concentrations of the solvent species ( S , Uf , and Vf) 
the dynamics approaches a reaction-diffusion model in which advection as well as momentum 
and heat transport become negligible. The reference (initial) values for mole fraction are 
used to prescribe Dirichlet boundary conditions for species on each side of the domain. This 
setup mimics open reservoirs in the form of permeable membranes [43]; note however that 
implementing boundaries that are also open for advective mass transport (i.e., inflow and 
outflow) is quite challenging [119] and not presently supported in our implementation. At 
the boundaries, the temperature is hxed at T = 300K and the fluid velocity is set to zero 
(no-slip) so species transport is primarily clue to mass diffusion. 

In our hydrodynamic simulations the fluid is modeled as a hard sphere dilute gas so the 
transport coefficients depend on the masses and diameters of the particles of each species. 
The particles for all species in the BPM model have equal mass (m = 6.64 • 1CT 26 g) so as to 
ensure that the reactions conserve mass. For all species the particles have only translational 
energy and no internal degrees of freedom (i.e., z = 0) so pressure and enthalpy are unaffected 
by reactions. In the BPM model species U plays the role of the “inhibitor” while species V 
is the “activator.” [112] Typically pattern formation occurs when the inhibitor diffuses faster 
than the activator so we set the diameter of species U particles to be smaller than that of V 
particles, specifically d\ = 0.125 nm and da = 0.5 nm. Since we take kf /cf (see Table III), 
the intermediary species, W, supports the activation of V and thus we set the diameters of 
V and W to be equal (this makes these two species hydrodynamically indistinguishable). 
The diameters of the other species (S, Uf , and Vf) are small, d 4 = d 5 = d e = 0.025 nm, 
so they diffuse rapidly from the boundaries and within the system. These specific values of 
the diameters are chosen such that the self-diffusion coefficient of S (and Uf, Vf) is roughly 
an order of magnitude larger than the diffusion coefficient of U, which is itself an order of 
magnitude larger than the diffusion coefficient of V (and W). 

The system is simulated in a rectangular domain that is divided into 128 x 128 x 1 cells 
with Ax = Ay = 100 nm. The magnitude of the noise is varied by varying the domain 
thickness, which was either Az = 100 nm (low noise) or 10 nm (high noise). The reference 
value for the number of molecules per cell for the species of interest, U, is O(10 4 ) in the 
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4.38- 10 18 


W 

3.8- 10 17 


5 

5.0 • 10 2 ° 


Uf 

1.33 • 10 2 ° 


Vf 

5.0 • 10 2 ° 


TABLE III. Tables of reference number densities (bottom) and reaction rate parameters (top) for 
BPM model used to create Fig. 6 , in CGS units. 


former case and O(10 3 ) in the latter. The total number of solvent molecules per cell is O(10 6 ) 
for weak noise and O(10 5 ) for strong noise. The time step is At = 10 ps, as determined 
from stability requirements for the explicit temporal integrator. 

Figure 6 illustrates the pattern formation observed in the system for low noise (top 
row), high noise (middle row), and deterministic evolution (bottom row) started from a 
perturbed initial condition generated by the high noise simulation (see figure caption). The 
boundaries take some time to influence the center of the domain, so in the center the 
reservoir species are depleted and the system moves toward chemical equilibrium. However 
the boundary continuously forces the system so eventually spotted patterns develop, starting 
near the boundary, eventually filling the system. The resulting patterns are qualitatively 
similar to the “A pattern” observed by Pearson [98] in the GS model. In simulations with 
only species diffusion (i.e., setting all other transport to zero) we find similar patterning, 
indicating that this system is well-approximated by reaction-diffusion due to very large 
solvent concentrations. 

In [110] Turing writes, “Another implicit assumption concerns random disturbing influ¬ 
ences. Strictly speaking one should consider such influences to be continuously at work. This 
would make the mathematical treatment considerably more difficult without substantially 
altering the conclusions.” However, we see from Fig. 6 that the evolution is qualitatively 
different when large spontaneous fluctuations are present (middle row), as compared to 
when they are absent (bottom row). Specifically, the speed at which the patterns develop 
and propagate is noticeably accelerated by the spontaneous fluctuations (top and middle 
row), though the patterns themselves are qualitatively unchanged in this particular case. 
Other studies using reaction-diffusion models and particle simulations have reached similar 
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FIG. 6. Time sequence of images of p\ (species U) at t = 10, 15, 20, 25, and 30 ps (time step 
At = 10 ps, for other parameters see Table III). The top row is a low noise case corresponding to 
a domain of thickness of 100 nm; the middle row is a high noise case corresponding to a domain 
of thickness 10 nm; and the bottom row is the deterministic evolution with an initial perturbation 
introduced by keeping high noise on up to time t = 2 ns, and then turning the noise off. The color 
scale spans from pi = 3.3 • 10 -4 (blue) to p\ = 2.0 • ICC 3 (red). 

conclusions [4-7]. In [120] the authors investigated the Gray-Scott model by RDME sim¬ 
ulations and concluded, ’Complex spatiotemporal patterns, including spiral waves, Turing 
patterns, self-replicating spots and others, which are not captured or correctly predicted by 
the deterministic reaction-diffusion equations, are induced by internal reaction fluctuations.’ 


V. CONCLUSIONS AND FUTURE WORK 

In this work we have formulated a fluctuating hydrodynamics model for chemically reac¬ 
tive ideal gas mixtures, and developed a numerical algorithm to solve the resulting system 
of stochastic partial differential equations. In our Langevin formalism, the stochastic mass, 
momentum and heat flux as well as the stochastic chemical production rate, are modeled 
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using uncorrelated white noise processes, and the local number densities are real variables. 
This is contrast to the more traditional chemical master equation (CME) description of re¬ 
actions that accounts for every individual reaction as a small jump of the (potentially very 
large) integer number of reactant molecules. We formulated the thermodynamic driving 
force for chemical reactions in agreement with nonlinear nonequilibrium thermodynamics 
[59] and considered two different Langevin formulations of the stochastic chemical produc¬ 
tion rate. The first formulation is based on the law of mass action cast in the GENERIC 
framework [62] and leads to a noise covariance that is the logarithmic mean (LME) of the 
forward and reverse production rates. This formulation is fully consistent with equilibrium 
statistical mechanics, more specifically, the resulting dynamics is time reversible (i.e., sat¬ 
isfies detailed balance) with respect to the Einstein distribution for a closed system. The 
second formulation is based on the chemical Langevin equation (CLE) [47, 57, 69], and while 
it is not consistent with equilibrium statistical mechanics this form has its own merits. 

We compared the two formulations on two chemical reaction networks for a well-mixed 
system, for both a simple dimerization reaction and a more complex network exhibiting 
bistability. We confirmed that at thermodynamic equilibrium the LME is more appropriate 
than the CLE, however, this is reversed for systems away from equilibrium, when compared 
with the CME. Not unexpectedly, neither is found to be entirely appropriate for describing 
rare events or large deviations from equilibrium. These examples remind us that a stochastic 
differential equation, which is a diffusion process, cannot be a uniformly accurate approxi¬ 
mation for the CME, which is a Poisson process; the large-deviation statistics for Poisson 
noise is different than that of Gaussian noise. To further complicate the picture, there are 
known examples in which the discrete (integer-valued) nature of molecular populations plays 
a key role, implying that descriptions using real-valued concentrations such as fluctuating 
hydrodynamics must fail. For example, in wave fronts of the Fisher, Kolmogorov, Petro- 
vski, Piskunov (FKPP) type it has been shown that the discreteness of the ME induces a 
logarithmic correction to the wave speed [121], similar to that observed when introducing a 
small cutoff in the leading edge of a FKPP front [122] 

Another alternative coarse-grained description of stochastic chemistry, which we did not 
consider in this work, is tau leaping [66, 123, 124], Tau leaping is usually seen as a numerical 
method for approximately solving the CME, and the CLE can be seen as an approximation 
to tau leaping in a specific (central) limit in which a Poisson and a Gaussian variable 
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become indistinguishable [66]; observe however that the two kinds of distributions always 
have different tails. A different and more interesting characterization of tau leaping is to 
see it instead as an alternative to the CLE that maintains the Poisson nature of the noise 
rather than replacing it with Gaussian noise. In the limit that the time step At —> 0, for 
Gaussian noise one gets the CLE, and for Poisson noise one gets, in principle, the CME. 
In this limit using the SSA algorithm is an efficient method to solve the CME exactly, 
and tau leaping is only useful as a numerical scheme for large At. As mentioned earlier, 
computational fluctuating hydrodynamics should only be considered useful (or even valid) 
when each update of the coarse-grained degrees of freedom involves an average over many 
molecular events, such as many molecular collisions for momentum and energy transport, or 
many reactive collisions for chemistry. In other words, to distinguish it from the CME and 
the associated SSA algorithm, for tau leaping one should choose the time step size in a way 
that ensures that many reactions occur in each reactive channel. In this sense, tau leaping 
can be seen as a coarse-graining in time of the CME jump process, and, when combined with 
spatial coarse-graining, has the potential to be a useful coarse-grained model that bypasses 
the need for Langevin models of chemistry in our numerical schemes for reactive fluctuating 
hydrodynamics. Additional studies are needed to access the accuracy of tau leaping in 
situations where Langevin descriptions do poorly and such investigations are in progress. 

As is well-known, the failure of Langevin approximations to describe large deviations is 
in fact closely connected to the fact that traditional linear nonequilibrium thermodynamics 
fails to describe chemical reactions because the entropy production rate is generally a non¬ 
quadratic function of the thermodynamic driving force (affinity). The mesoscopic Kramers 
picture of chemical reactions, as developed for isomerization in Ref. [59], is an interesting 
approach, which, however, remains mostly of theoretical utility; numerical simulations of 
this model would need to handle additional dimensions, as well as very slow diffusion across 
the reaction barrier. Furthermore, it is not obvious how to extend this description to general 
multispecies fluids with complex reaction networks. 

While the CME is a well-agreed upon and well-justified description of statistically ho¬ 
mogeneous systems, the story is much less clear for systems with spatial inhomogeneity; 
in fact, the precise mathematical meaning of “well-mixed” and the range of validity of the 
RDME remains obscure. [21-23] A law of large numbers has been rigorously proven for sev¬ 
eral particle models and takes the expected form of a deterministic reaction-diffusion partial 
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differential equation. [72, 73] Regarding fluctuations, however, there are very few particle 
models for which central limit theorems [125] or large deviations functionals [126] are known 
explicitly. It has been demonstrated that inhomogeneity leads to qualitative changes in the 
nature of phase transitions in bistable systems [127]. It is also known that fluctuations can 
effectively renormalize the macroscopic transport and lead to non-analytic corrections to 
the law of mass action [2, 3]. It remains to be seen whether nonlinear spatially-extended 
fluctuating hydrodynamics models can correctly reproduce this effect compared to particle 
simulations. 15 

Furthermore, most particle models used for reaction-diffusion problems assume indepen¬ 
dent random Brownian walkers that react when coming near each other, and thus completely 
neglect transport mechanisms such as advection, sound waves, thermal conduction, etc. Fur¬ 
thermore, the mechanism of diffusion used in these models implicitly neglects the long-ranged 
hydrodynamic correlations present among particles diffusing in a viscous solvent [52], No¬ 
table exceptions are variants of the Direct Simulation Monte Carlo method (DSMC), which 
use a dilute [128, 129] or dense [130] gas kinetic theory description of momentum and energy 
transport fully consistent with fluctuating hydrodynamics [130]. While chemical reactions 
are commonly included in DSMC schemes [131, 132], further studies of fluctuations and their 
consistency with nonequilibrium thermodynamics are needed. While some investigations of 
spatially-distributed reactive systems have been performed using DSMC [118], a careful 
comparison to coarse-grained mesoscopic descriptions such as fluctuating hydrodynamics is 
needed. To model realistic chemistry modern DSMC codes use either the Total Collision 
Energy (TCE) model [128] or the more recent Quantum-Kinetic (QK) model [133], both of 
which we plan to compare with our formulations of reactive fluctuating hydrodynamics. 

In Section IV B we studied the coupling of velocity fluctuations and chemistry in a sys¬ 
tem kept in a non-equilibrium steady state via boundary conditions. We found that, in 
agreement with existing theoretical computations, the chemical reactions have a strong ef¬ 
fect on the giant long-range correlated concentration fluctuations. In this work we focused 
on gas mixtures, for which the Schmidt number is of order unity. We found that reactions 
profoundly change the nature of the giant fluctuations, whose spectrum switches from the 
well-known ~ k~ A at large wavenumbers to a flat plateau (with a value controlled by the 
15 Fluctuating hydrodynamics simulations have been shown to correctly describe the renormalization of 
diffusion coefficients in binary fluid mixtures [105], but we are not aware of any work in this direction for 
reactive systems. 
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reaction rate) at small wavenumbers. This could be useful, for example, for experimentally 
measuring reaction rates. However, we found that the simple quasi-periodic near-equilibrium 
theory we constructed was not in quantitative agreement with the simulations, indicating 
that a more precise theory is needed. 

Reactive liquid mixtures are common in practice and exhibit interesting coupling of hydro¬ 
dynamics with transport that has been studied both theoretically [134] and experimentally 
[9]. The Schmidt number in liquids is, however, very large, and the compressibility is very 
small (i.e., the speed of sound is very large), making the method presented here compu¬ 
tationally infeasible. In the future, we will consider extending low Mach number (quasi- 
incompressible) methods that treat momentum diffusion implicitly to include stochastic 
chemistry, by combining Langevin or tail-leaping based descriptions of chemistry with the 
formulation and numerical methods developed in Refs. [44, 135] for general non-ideal liquid 
mixtures. 

Finally, the influence of hydrodynamic fluctuations on reactions will likely be very im¬ 
portant for surface chemistry. [136, 137] An important application is heterogeneous catalysis 
in which a highly reactive catalytic surface facilitates bond breaking and bond rearrange¬ 
ment of adsorbed molecules. In this context mesoscale simulations are particularly useful 
for the study of nano-catalytic systems [138]. Microscopic catalytic particles have many 
advantages, such as higher activity, increased selectivity, and longer lifetime. However, to 
operate effectively these particles must have electrical contact with a substrate (typically a 
flat surface) and they must not be so small as to not have enough electrons to catalyze a 
reaction. For example, nanowire catalysts are typically 10-100 nm in size so the hydrody¬ 
namic environment in which the chemistry and transport occur is of mesoscopic scale (a few 
microns). 

As in the case with chemistry in bulk flow, particle-based simulations will be useful bench¬ 
marks for comparison with reactive fluctuating hydrodynamics modeling surface chemistry. 
Furthermore, molecular simulations can be embedded within a fluctuating hydrodynamic 
code to create an Algorithm Refinement (AR) hybrid. [39, 139] The idea is to use a particle- 
based simulation for the domain near the surface to capture the physics at the molecular 
scale, such as adsorption of reactants onto the surface, surface diffusion, surface reactions, 
and desorption of products. This particle-based simulation would then be coupled to a fluc¬ 
tuating hydrodynamic simulation that treats the bulk fluid. Previous work for non-reactive 
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fluids has already proven the utility of AR hybrids for simple interfaces [53] and this numer¬ 
ical framework should prove useful in the study of surface reactions and active membrane 
transport. 
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Appendix A: Deterministic and Stochastic Transport in Ideal Gas Mixtures 

This appendix summarizes fluctuating hydrodynamics for ideal gas mixtures; for a more 
general and detailed exposition see [42], Each of the hydrodynamic transport terms in (1)- 
(3) contains a deterministic term, denoted with an overbar and a stochastic term denoted 
by a tilde (e.g., T = T + T). The deterministic viscous tensor is, 

where r\ and n are the shear and bulk viscosity, respectively. We neglect any possible effect of 
the chemical reactions on the transport coefficients of the mixture. For example, we neglect 
any coupling between bulk viscosity and chemical reactions, which, in principle is allowed by 
the Curie principle since both are scalar processes. [74] The corresponding stochastic viscous 
flux tensor can be written as [36, 140] 

n(r,i)= A^Z n +(' V ^- V /^i')T 1 (Z"). (A2) 

The symmetric Gaussian random tensor field, Z n , is formulated as Z u = (Z n + ( Z n ) T ') /y/2 
where Z n is a white-noise random Gaussian tensor field, that is, (Z n ) = 0 and, 

(^5( r > t ) Z 5( r/ ’ 0) = 5 onfysS{r ~ r') 5(t - t') . 


n = -7] (Vv + (Vv) T ) 


g 7 ]I(V-v), 


(Al) 
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with a, (3, 7 , 8 — {x, y , z} being spatial components. 

The deterministic mass flux and heat flux depend on the gradients of concentration, 
pressure, and temperature. For an ideal gas we can write the deterministic fluxes as 

XT' 


t = - P yr> (d + xx 


T 


and 


(A3) 

(A4) 


Q = -AVT + (ksTx 1 AT 1 + h 1 )T 

where h is a vector of specific enthalpies. Here the diffusion driving force [83] is 

d = VX + (X - 1") —. 

P 

and X, y and JY[ are diagonal matrices of mole fractions X, mass fractions Y and molecular 
masses M. The matrix of multicomponent flux diffusion coefficients, D, the vector of rescaled 
thermal diffusion ratios, x, and the thermal conductivity, A, can be obtained from standard 
software libraries, such as EGLIB [141], or from standard references, such as Hirshfelder et 
al. [142], 

For the ideal gas transport coefficients described above, we define 


L = p yyoy, 
ks 


£ = ksTM-'x, C = T z A. 


(AS) 


where m = Y s /m s )~ l is the mixture-averaged molecular mass. The stochastic terms for 
the combined species equations and energy equation are determined from the phenomeno¬ 
logical equations of nonequilibrium thermodynamics that relate fluxes to thermodynamic 
driving forces through the Onsager matrix, and the fluctuation-dissipation balance princi¬ 
ple. [42] Specifically, the stochastic mass flux is 

F = B Z T (A6) 

where BB 7 = 2/q,L, and Z T is a white-noise random Gaussian vector field with uncorrelated 
components, that is, (. Z T ) = 0 and 

( z as( r ^) z y s '( r',0) = <*(r “ r ') S(t - t') 

with a, (3 = {a;, y, z} being spatial components. Note that the matrix B can be obtained 
using Cholesky decomposition; the constraint T s = 0 is ensured by construction. [42] 
The stochastic heat flux is 

Q = + h T )F, (A7) 

where (Z®(r, t)Zf(r', t')) = 5 aP 5(r - r') S(t - t'). 
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Appendix B: Equilibrium distribution for a dimerization reaction 


The entropy of mixing for a system undergoing dimerization is S — k B hi N p , where N p 
is the number of distinct ways of forming M dimers out of a total of N monomers. This is 
straightforward to compute. The number of ways to choose 2M out of the N atoms to be 


in pairs is, 

/ N \ _ N\ 

\2M ) ~ (2 M)\(N - 2 M)\ 

Next we need to group the 2 M atoms into pairs; the number of distinct ways of pairing 2 M 
objects is (2M)\/2 m M\. This gives the entropy 

S = kB hl ( 2 M M! (iV~— 2M)l ) + k ^ M ’ (B1) 

where we have included an additional reference chemical potential /i to set the equilibrium 
concentration. 

By expanding the logarithm of the right-hand side of (Bl) using Stirling’s leading-order 
approximation we obtain the thermodynamic limit of the entropy as a function of the 
monomer mass fraction Y\ = (. N — 2M) /N. After fixing the chemical potential from the 
requirement that the most probable mass fraction is (Yi) eq = 1/2, we get, 

S = Nk B ~ In (1 - Y,) (1 - Tj) - Y 1 In (Tj) - i Tj (In (2) - 1) , (B2) 

This gives an Einstein distribution P ~ e s / kB exactly matching (27). 

It is useful to compare this equilibrium distribution of the LME with that of the CME for 
a well-mixed system of volume AV. It is not hard to show that the CME for a dimerization 
reaction is in detailed balance with respect to the Einstein distribution with entropy (Bl), 
with the reference chemical potential set to 

" = 111 (^7f) ' (B3) 

We can use this to construct a rather accurate continuum approximation to this exact 


microscopic result by including the next order term in the Stirling formula, 


N\ « V2tcN 



(B4) 


when expanding (Bl). This gives a better continuous approximation (labeled “Stirling” in 
the figures in the main body of this paper) to the discrete Einstein distribution than (B2), 
especially for small number of particles in the well-mixed cell. 
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Appendix C: Giant fluctuations in the presence of reactions 


This appendix outlines the fluctuating hydrodynamics theory for the long-range correla¬ 
tions of concentration fluctuations in a binary mixture undergoing a dimerization reaction, 
as studied numerically in Section IVB. We neglect the Dufour effect and assume the system 
to be isothermal, taking contributions from temperature fluctuations to be of higher order. 
Furthermore, we neglect gravity, assume the system is incompressible, and take the density 
and transport coefficients to be constant. We consider a “bulk” system [35], i.e., we neglect 
the influence of the boundaries. This gives an accurate approximation for wavenumbers 
that are large compared to the inverse height of the domain; for smaller wavenumbers the 
boundaries are expected to suppress the giant fluctuations [35, 109]. We also neglect the 
stochastic mass flux in the concentration equation since we are concerned with the nonequi¬ 
librium contribution due to the forcing by the velocity fluctuations. 

We assume that all of the concentration gradients are in the same direction (say, the 
y axis). The incompressibility constraint is most easily handled by applying a V x Vx 
operator to the momentum equation [35] to obtain a system involving only the component of 
the velocity parallel to the gradient, vy = v y . The same calculation can easily be generalized 
to a multispecies mixture as well, see Appendix B in Ref. [42], This system of equations 
can be most easily solved in the Fourier domain, by considering wavevectors k in the plane 
perpendicular to the gradient, k = k±. 

It is very straightforward to derive (29) in the limit of large Schmidt number by con¬ 
sidering the fluctuating concentration equation forced by an overdamped (steady Stokes) 
fluctuating velocity. The steady Stokes equation in Fourier space has the form, 


r)k 2 v\\ = s/2)]k B T 0 ikW(t), 


where To is the constant temperature. Here W(f) is a white-noise process (one per wavenum¬ 
ber), giving the white-in-time velocity v\\(t) = k~ l sJ^kjfTfffr] iW{t). A general form of the 
linearized equation for the concentration fluctuations, c = 8Y\ = Y\ — (Yi), in Fourier space 
will have the form 


d t c = — Dk 2 c — f>c — v\\f = — (Dk 2 + ip) c — if y 


2k B T 0 
r) k 2 


Wit), 


where is the reaction rate at the equilibrium point (around which we are linearizing), and 
f = d (I'i) /dy is the applied concentration gradient. For our specific reaction and the choice 


44 



of equilibrium point (Yi) eq = 1/2, we have that ^ = 3 k~ (to see this simply linearize (31) 
around {Y\ ) eq = 1/2). The resulting nonequilibrium static structure factor (recall that here 
we neglect the contribution due to the stochastic mass flux) is straightforward to calculate 
(see, for example, Appendix B in Ref. [42]), 


S neq (k) = <cc*> = 




-f 2 = 




7]k 2 ( Dk 2 + ^) J r]Dk A (1 + ipD^k- 2 ) 
which is exactly (29) with the penetration depth 




(Cl) 


d 2 = 


D 


D 

3jF' 


Since the Schmidt number is not very large for gases, we should improve the theory by 
not taking the overdamped limit but rather adding a velocity equation and considering the 
linearized inertial equations, 


d t c = —Dk 2 c — i\)c — v\\f 
Pod t v\\ = -r]k 2 v || + ^/2r]k B T 0 ikW(t), 


where po is the equilibrium density and v = rj/po is the kinematic viscosity. Solving this 
system of linear SODEs gives the improved structure factor 


Sneq(k) 


ksTo 


r]Dk A (1 + ^D~ l k~ 2 ) + ^v^k- 2 ) 


f 


(C2) 


which shows that the finite value of the Schmidt number S c = v/D has an important effect 
on the giant fluctuations. In particular, in the more complete theory (C2) 


s -( fc = °> = ^ /2 

is finite and not infinite as in (Cl), which assumes infinite Schmidt number. 
In the more complete theory we see the following three regimes: 


1. For large wavenumbers we have 


S', 


neq 



ksTp 2 
r]Dk AJ 


as if there were no reaction. 
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2. If the Schmidt number is small, D ~ u, as for gases, then for small wavenumbers we 
instantly switch to a flat spectrum 


5 , 


neq 




fcgTo ,2 


3. If the Schmidt number S'c is large, D <C zy as for liquids, then for intermediate 
wavenumbers 



ksTp ,2 
ij'ipk 2 


we observe a change in the power law to S neq (k) ~ k~' 2 . Note however that this range 
spans only \/S c orders of magnitude in k, so even for S c ~ 10 4 (which applies to 
macromolecular solutions), the k~ 2 power-law only extends over at most two decades. 
Therefore, even for liquids with large Schmidt numbers the inertial equations should 
be used to model giant fluctuations in reactive mixtures. 


[1] A. Lemarchand and B. Nowakowski, Molecular Simulation 30, 773 (2004). 

[2] K. Kang and S. Redner, Phys. Rev. Lett. 52, 955 (1984). 

[3] A. A. Winkler and E. Frey, Phys. Rev. Lett. 108, 108301 (2012). 

[4] A. Lemarchand and B. Nowakowski, EPL (Europhysics Letters) 94, 48004 (2011). 

[5] P. Dziekan, L. Signon, B. Nowakowski, and A. Lemarchand, The Journal of chemical physics 
139, 114107 (2013). 

[6] Y. Cao and R. Erban, Bulletin of mathematical biology 76, 3051 (2014). 

[7] F. Zheng-Ping, X. Xin-Hang, W. Hong-Li, and O. Qi, Chinese Physics Letters 25, 1220 
(2008). 

[8] P. Dziekan, J. S. Hansen, and B. Nowakowski, The Journal of Chemical Physics 141, 124106 
(2014). 

[9] C. Almarcha, P. M. Trevelyan, P. Grosfils, and A. De Wit, Physical Review E 88, 033009 
(2013). 

[10] D. Brogioli, Phys. Rev. Lett. 105, 058102 (2010). 

[11] D. Brogioli, Phys. Rev. E 84, 031931 (2011). 

[12] D. Brogioli, The Journal of Chemical Physics 139, 184102 (2013). 


46 



[13] A. Mahmutovic, D. Fange, O. G. Berg, and J. Elf, Nature methods 9, 11631166 (2012). 

[14] M. Malek-Mansour and G. Nicolis, Journal of Statistical Physics 13, 197 (1975). 

[15] G. Nicolis and I. Prigogine, Self-organization in nonequilibrium systems: from dissipative 
structures to order through fluctuations, A Wiley-Interscience Publication (Wiley, 1977). 

[16] C. W. Gardiner, Handbook of stochastic methods: for physics, chemistry & the natural sci¬ 
ences, 3rd ed., Series in synergetics, Vol. Vol. 13 (Springer, 2003). 

[17] R. Erban, J. Chapman, and P. Maini, arXiv preprint arXiv:0704.1908 (2007). 

[18] D. Fange, A. Mahmutovic, and J. Elf, Bioinformatics 28, 3155 (2012), 

http://bioinformatics, oxfordjournals. org/content/28/23/3155. full, pdf+html. 

[19] D. T. Gillespie, Journal of Computational Physics 22, 403 (1976). 

[20] D. T. Gillespie, Annual Review of Physical Chemistry 58, 35 (2007). 

[21] D. Fange, O. G. Berg, P. Sjoberg, and J. Elf, Proceedings of the National Academy of 
Sciences 107, 19820 (2010). 

[22] S. Hellander, A. Hellander, and L. Petzold, Phys. Rev. E 91, 023312 (2015). 

[23] S. A. Isaacson, The Journal of Chemical Physics 139, 054101 (2013). 

[24] M. Dobrzynski, J. V. Rodrguez, J. A. Kaandorp, and J. G. Blorn, Bioinformatics 23, 1969 
(2007), http://bioinformatics.oxfordjournals.org/content/23/15/1969.full.pdf+html. 

[25] S. S. Andrews, N. J. Addy, R. Brent, and A. P. Arkin, PLoS Comput Biol 6, el000705 

( 2010 ). 

[26] J. S. van Zon and P. R. ten Wolde, Phys. Rev. Lett. 94, 128103 (2005). 

[27] J. S. van Zon and P. R. ten Wolde, J. Chem. Phys. 123, 234910 (2005). 

[28] T. Oppelstrup, V. V. Bulatov, A. Donev, M. H. Kalos, G. H. Gilmer, and B. Sadigh, Phys. 
Rev. E 80, 066701 (2009). 

[29] A. Donev, V. V. Bulatov, T. Oppelstrup, G. H. Gilmer, B. Sadigh, and M. H. Kalos, J. 
Comp. Phys. 229, 3214 (2010). 

[30] A. J. Mauro, J. K. Sigurdsson, J. Shrake, P. J. Atzberger, and S. A. Isaacson, J. Comp. 
Phys. 259, 536 (2014). 

[31] A. Bezzola, B. B. Bales, R. C. Alkire, and L. R. Petzold, J. Comp. Phys. 256, 183 (2014). 

[32] D. T. Gillespie, E. Seitaridou, and C. A. Gillespie, The Journal of Chemical Physics 141, 
234115 (2014). 

[33] K. Kuo, Principle of Combustion (John Wiley & Sons, 2005). 


47 



[34] C. Law, Combustion Physics (Cambridge University Press, 2006). 

[35] J. M. O. D. Zarate and J. V. Sengers, Hydrodynamic fluctuations in fluids and fluid mixtures 
(Elsevier Science Ltd, 2007). 

[36] L. Landau and E. Lifshitz, Fluid Mechanics , Course of Theoretical Physics, Vol. 6 (Pergamon 
Press, Oxford, England, 1959). 

[37] A. Garcia, M. Malek-Mansour, G. Lie, M. Mareschal, and E. Clementi, Phys. Rev. A 36, 
4348 (1987). 

[38] R. Adhikari, K. Stratford, M. E. Cates, and A. J. Wagner, Europhysics Letters 71, 473 
(2005), arXiv:cond-mat/0402598. 

[39] J. Bell, A. Garcia, and S. Williams, ESAIM: M2AN 44, 1085 (2010). 

[40] A. Donev, E. Vanden-Eijnden, A. L. Garcia, and J. B. Bell, Communications in Applied 
Mathematics and Computational Science 5, 149 (2010). 

[41] B. Shang, N. Voulgarakis, and J. Chu, J. Chem. Phys. 137, 044117 (2012). 

[42] K. Balakrishnan, A. L. Garcia, A. Donev, and J. B. Bell, Phys. Rev. E 89, 013017 (2014). 

[43] A. Donev, A. J. Nonaka, Y. Sun, T. G. Fai, A. L. Garcia, and J. B. Bell, Communications 
in Applied Mathematics and Computational Science 9, 47 (2014). 

[44] A. Donev, A. J. Nonaka, A. K. Bhattacharjee, A. L. Garcia, and J. B. Bell, “Low Mach 
Number Fluctuating Hydrodynamics of Multispecies Liquid Mixtures,” (2015), to appear in 
Physics of Fluids, ArXiv:1412.6503. 

[45] W. Koh and K. T. Blackwell, The Journal of Chemical Physics 137, 154111 (2012), 
http://dx.doi.Org/10.1063/l.4758459. 

[46] A. Ghosh, A. Leier, and T. T. Marquez-Lago, Theoretical Biology and Medical Modelling 
12, 5 (2015). 

[47] D. T. Gillespie, The Journal of Chemical Physics 113, 297 (2000). 

[48] P. N. Segre, R. W. Gammon, and J. V. Sengers, Phys. Rev. E 47, 1026 (1993). 

[49] A. Vailati and M. Giglio, Nature 390, 262 (1997). 

[50] A. Vailati, R. Cerbino, S. Mazzoni, C. J. Takacs, D. S. Cannell, and M. Giglio, Nature 
Communications 2, 290 (2011). 

[51] A. Donev, T. G. Fai, and E. Vanden-Eijnden, Journal of Statistical Mechanics: Theory and 
Experiment 2014, P04004 (2014). 

[52] A. Donev and E. Vanden-Eijnden, J. Chem. Phys. 140, 234115 (2014), 


48 



http://dx.doi.org/10.1063/L4883520. 

[53] A. Donev, J. B. Bell, A. L. Garcia, and B. J. Alder, SIAM J. Multiscale Modeling and 
Simulation 8, 871 (2010). 

[54] B. Franz, M. B. Flegg, S. J. Chapman, and R. Erban, SIAM Journal on Applied Mathematics 
73, 1224 (2013). 

[55] H. Haken, Synergetics: Introduction and Advanced Topics, Physics and Astronomy Online 
Library (Springer, 2004). 

[56] P. J. Atzberger, J. Comp. Phys. 229, 3474 (2010). 

[57] J. Keizer, Statistical thermodynamics of nonequilibrium processes (Springer, 1987). 

[58] I. Pagonabarraga, A. Perez-Madrid, and J. Rubi, Physica A: Statistical Mechanics and its 
Applications 237, 205 (1997). 

[59] D. Bedeaux, I. Pagonabarraga, J. de Zarate, J. Sengers, and S. Kjelstrup, Phys. Chern. 
Chem. Phys. 12, 12780 (2010). 

[60] J. de Zarate, J. Sengers, D. Bedeaux, and S. Kjelstrup, J. Chem. Phys. 127, 034501 (2007). 

[61] D. Bedeaux, J. O. de Zarate, I. Pagonabarraga, J. Sengers, and S. Kjelstrup, J. Chem. Phys. 
135, 124516 (2011). 

[62] H. Ottinger and M. Grmela, Physical Review E 56, 6633 (1997). 

[63] P. Espanol, J. Anero, and I. Zuniga, J. Chem. Phys. 131, 244117 (2009). 

[64] J. de la Torre, P. Espaol, and A. Donev, 142, 094115 (2015), 

http://dx.doi.Org/10.1063/l.4913746. 

[65] R. Grima, P. Thomas, and A. V. Straube, The Journal of chemical physics 135, 084103 

( 2011 ). 

[66] D. Gillespie, A. Hellander, and L. Petzold, J. Chem. Phys. 138 (2013). 

[67] H. Grabert, P. Hanggi, and I. Oppenheim, Physica A: Statistical Mechanics and its Appli¬ 
cations 117, 300 (1983). 

[68] P. Hanggi, H. Grabert, P. Talkner, and H. Thomas, Phys. Rev. A 29, 371 (1984). 

[69] D. F. Anderson and T. G. Kurtz, in Design and analysis of biomolecular circuits (Springer, 
2011) pp. 3-42. 

[70] F. Baras, J. E. Pearson, and M. M. Mansour, The Journal of Chemical Physics 93, 5747 
(1990). 

[71] F. Baras, M. M. Mansour, and J. E. Pearson, The Journal of Chemical Physics 105, 8257 


49 



(1996). 

[72] L. Arnold and M. Theodosopulu, Advances in Applied Probability 12, pp. 367. 

[73] L. Arnold, in Dynamics of synergetic systems (Springer, 1980) pp. 107-118. 

[74] S. Grossmann, The Journal of Chemical Physics 65, 2007 (1976). 

[75] H. C. Ottinger, (2015), draft. 

[76] D. Bedeaux, S. Kjelstrup, and H. C. Ottinger, The Journal of chemical physics 141, 124102 
(2014). 

[77] D. Bedeaux and S. Kjelstrup, Physical Chemistry Chemical Physics 10, 7304 (2008). 

[78] H. C. Ottinger, Beyond equilibrium thermodynamics (Wiley Online Library, 2005). 

[79] Y. L. Klimontovich, Physics-Uspekhi 37, 737 (1994). 

[80] M. Hritter and H. Ottinger, J. Chern. Soc., Faraday Trans. 94, 1403 (1998). 

[81] H. Grabert, Projection operator techniques in nonequilibrium statistical mechanics (Springer- 
Verlag, 1982). 

[82] S. Kjelstrup, D. Bedeaux, E. Johannessen, and J. Gross, Non-Equilibrium Thermodynamics 
for Engineers (World Scientific Pub Co Inc, 2010). 

[83] G. D. C. Kuiken, Thermodynamics of Irreversible Processes: Applications to Diffusion and 
Rheology (Wiley, 1994). 

[84] D. Kondepudi and I. Prigogine, Modern thermodynamics: from heat engines to dissipative 
structures (John Wiley, 1998). 

[85] M. Zemansky and R. Dittman, Heat and Thermodynamics: An Intermediate Textbook , In¬ 
ternational student edition: MacGraw-Hill (McGraw-Hill, 1981). 

[86] R. Pathria and P. Beale, Statistical Mechanics (Elsevier Science, 2011). 

[87] V. Giovangigli, Multicomponent Flow Modeling (Birkhauser Boston, 1999). 

[88] A. Mielke, Calculus of Variations and Partial Differential Equations 48, 1 (2013). 

[89] A. Mielke, Nonlinearity 24, 1329 (2011). 

[90] S. Delong, B. E. Griffith, E. Vanden-Eijnden, and A. Donev, Phys. Rev. E 87, 033302 (2013). 

[91] S. Delong, Y. Sun, , B. E. Griffith, E. Vanden-Eijnden, and A. Donev, Phys. Rev. E 90, 
063312 (2014). 

[92] F. B. Usabiaga, J. B. Bell, R. Delgado-Buscalioni, A. Donev, T. G. Fai, B. E. Griffith, and 
C. S. Peskin, SIAM J. Multiscale Modeling and Simulation 10, 1369 (2012). 

[93] P. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Stochastic 


50 



Modelling and Applied Probability (Springer Berlin Heidelberg, 2010). 

[94] H. Ottinger, Stochastic Processes in Polymeric Fluids: Tools and Examples for Developing 
Simulation Algorithms (Springer Berlin Heidelberg, 2012). 

[95] D. Anderson and J. Mattingly, Communications in Mathematical Sciences 9, 301 (2011). 

[96] P. Gray and S. Scott, Chemical Engineering Science 38, 29 (1983). 

[97] P. Gray and S. K. Scott, The Journal of Physical Chemistry 89, 22 (1985). 

[98] J. E. Pearson, 261, 189 (1993). 

[99] N. G. van Kampen, Stochastic Processes in Physics and Chemistry , 3rd ed. (Elsevier, 2007). 

[100] S. N. Ethier and T. G. Kurtz, Markov processes: characterization and convergence , Vol. 282 
(John Wiley & Sons, 2009). 

[101] D. T. Gillespie, The Journal of Physical Chemistry A 106, 5063 (2002). 

[102] A. Duncan, S. Liao, T. Vejchodsky, R. Erban, and R. Grirna, arXiv preprint arXiv: 1407.8256 
(2014). 

[103] M. Heymann and E. Vanden-Eijnden, Communications on pure and applied mathematics 
61, 1052 (2008). 

[104] A. Donev, A. L. Garcia, A. de la Fuente, and J. B. Bell, Phys. Rev. Lett. 106, 204501 

( 2011 ). 

[105] A. Donev, A. L. Garcia, A. de la Fuente, and J. B. Bell, J. of Statistical Mechanics: Theory 
and Experiment 2011, P06014 (2011). 

[106] H. N. W. Lekkerkerker and W. G. Laidlaw, Phys. Rev. A 9, 346 (1974). 

[107] J. L. Hita and J. M. O. de Zarate, Physical Review E 87, 052802 (2013). 

[108] S. R. DeGroot and P. Mazur, Non-Equilibrium Thermodynamics (North-Holland Publishing 
Company, Amsterdam, 1963). 

[109] J. M. Ortiz de Zarate, J. A. Fornes, and J. V. Sengers, Phys. Rev. E 74, 046305 (2006). 

[110] A. M. Turing, Philosophical Transactions of the Royal Society of London B: Biological Sci¬ 
ences 237, 37 (1952). 

[111] S. Kondo and T. Miura, 329, 1616 (2010). 

[112] M. Cross and H. Greenside, Pattern Formation and Dynamics in Nonequilibrium Systems 
(Cambridge University Press, 2009). 

[113] M. Malek-Mansour and J. Houard, Physics Letters A 70, 366 (1979). 

[114] T. Biancalani, D. Fanelli, and F. Di Patti, Phys. Rev. E 81, 046215 (2010). 


51 



[115] D. Fanelli, C. Cianci, and F. Di Patti, The European Physical Journal B 86, 1 (2013). 

[116] N. Kumar and W. Horsthemke, Phys. Rev. E 83, 036105 (2011). 

[117] E. P. Zemskov, K. Kassner, M. J. B. Hauser, and W. Horsthemke, Phys. Rev. E 87, 032906 
(2013). 

[118] P. Dziekan, A. Lemarchand, and B. Nowakowski, The Journal of Chemical Physics 137, 
074107 (2012). 

[119] R. Delgado-Buscalioni and A. Dejoan, Phys. Rev. E 78, 046708 (2008). 

[120] H. Wang, Z. Fu, X. Xu, and Q. Ouyang, The Journal of Physical Chemistry A 111, 1265 
(2007). 

[121] J. S. Hansen, B. Nowakowski, and A. Lemarchand, The Journal of Chemical Physics 124, 
034503 (2006). 

[122] E. Brunet and B. Derrida, Phys. Rev. E 56, 2597 (1997). 

[123] Y. Hu, T. Li, and B. Min, J. Chern. Phys. 135, 024113 (2011). 

[124] D. F. Anderson and M. Koyarna, Multiscale Modeling & Simulation 10, 1493 (2012). 

[125] C. Boldrighini, A. De Masi, and A. Pellegrinotti, Stochastic processes and their applications 
42, 1 (1992). 

[126] G. Jona-Lasinio, C. Landirn, and M. Vares, Probability theory and related fields 97, 339 
(1993). 

[127] S. Tanase-Nicola and D. K. Lubensky, Physical Review E 86, 040103 (2012). 

[128] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Clarendon, Ox¬ 
ford, 1994). 

[129] F. J. Alexander and A. L. Garcia, Computers in Physics 11, 588 (1997). 

[130] A. Donev, A. L. Garcia, and B. J. Alder, Phys. Rev. Lett 101, 075902 (2008). 

[131] M. J. Goldsworthy, M. N. Macrossan, and M. M. Abdel-Jawad, Physics of Fluids 19, 066101 
(2007). 

[132] C. Lilley and M. Macrossan, Physics of Fluids 16, 2054 (2004). 

[133] G. A. Bird, Physics of Fluids (1994-present) 23, 106101 (2011). 

[134] P. M. J. Trevelyan, C. Almarcha, and A. De Wit, Phys. Rev. E 91, 023001 (2015). 

[135] A. J. Nonaka, Y. Sun, J. B. Bell, and A. Donev, “Low Mach Number Fluctuating Hy¬ 
drodynamics of Binary Liquid Mixtures,” (2015), submitted to CAMCOS, ArXiv preprint 
1410.2300. 


52 



[136] G. Somorjai and Y. Li, Introduction to Surface Chemistry and Catalysis (John Wiley & Sons, 

2010 ). 

[137] G. A. Somorjai and Y. Li, Proceedings of the National Academy of Sciences 108, 917 (2011). 

[138] A. Z. Moshfegh, Journal of Physics D: Applied Physics 42, 233001 (2009). 

[139] S. A. Williams, J. B. Bell, and A. L. Garcia, SIAM J. Multiscale Modeling and Simulation 
6, 1256 (2008). 

[140] P. Espanol, Physica A 248, 77 (1998). 

[141] A. Ern and V. Giovangigli, EGLIB: A general-purpose Fortran library for multicomponent 
transport property evaluation, Version 3-4 (2004). 

[142] J. Hirschfelder, C. Curtiss, and R. Bird, Molecular Theory of Gases and Liquids, Structure 
of matter series (Wiley, 1967). 


53 



